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All About Fourier Analysis 

The first installment in a new series exploring the worid of time, 
firequency, and waveforms. 



Is there an engineer alive who 
hasn't heard the terms Fourier 
series, Fotirier transform, and 
Fovirier analysis? Probably not. 
Even if you haven't used the 
techniques in quite a while, if ever, yoa 
probably remember that it has some- 
thing to do with sine waves. But that 
simple statement seems entirely too 
innocuous to indicate the profound 
effect that Fourier analysis has had on 
almost every aspect of raigineering, 
and on everyday life. 

In this two-part series, we'll be look- 
ing at die motivation and concq)ts 
behind Fourier analysis, and why it's 
so important. This month, we'll look at 
the Fourier series, which allows us to 
represent a periodic function by a sum 
of many different sine waves. In the 
process, I'll give you a neat and sim- 
ple, but powerful, little algorithm of 
mine, which I call the "half-fast 
Fourier transform." Next time, we'll 
extend these concepts of sums and har- 
monics to the realm of infinities and 
infinitesimals, which will lead us to the 
all-important Fourier transform and its 
implementation in the famous Fast 
Fourier l^aasfcmn (FFT). 

getung motivated 

When I was in school, I found 
that I always got more out 
of a class if I knew what 
direction the professor's lectures were 
taking, and why. Whenever I set out to 
explain a new technique, I try to pro- 
vide the reasoning behind it, so readers 
can understand the background and 
usefulness of the technique. Basically, 



might have been going thi?mig^ die 
head of the inventor. 

Providing the reasoning for some 
topics is easier than others. Few have 
such a wealth of background as Fourier 
analysis; in Has case, we have an 
embarrassment of riches. Though you 
may not realize it, the concepts we 
need for the subject of Fourier an^ysis 
are already as familiar to you as your 
thumb and forefinger. 

When you were shcqiping for your 
first stereo system, you no doubt 
already knew that if an audio system is 
to faithfully reproduce music or other 
sounds, it must have "flat" frequency 
response. That is, it must amplify all 
frequencies in the audio range (20 Hz 
to 20,000 Hz) with equal amplitude. If 
it doesn't, the music won't soimd right. 
You also knew fliat every amplifier has 
a pair of tone control knobs labeled 
"bass" and "treble." Turning the bass 
up or down affects the reproduction of 
low frequencies, while the treble knob 
affects high frequencies. Turn the tre- 
ble up and the bass down, and the 
music sounds thin and sibilant. Turn 
them the opposite way and the music 
soimds muddy and muffled with 
booming bass drums. 

When you first tweaked one of those 
knobs yourself, you may have been 
siuprised to discover that these knobs 
made a profound difference, even if 
only one instrument, say a violin, were 
playing. That's because all musical 
instruments generate overtones, or har- 
monics, which are sounds with fre- 
quencies above the fundamental tone 
that defines the note being played. 
I^djusttng the tone contiols affecte the 





Matiiematically, 
the Fourier series 
is oniy a subset of 
a iarger set of 
metliods based on 



functions. 



relative amplitudes between the funda- 
mental and its overtones, which 
changes the sound quality. 

Implicit in all this talk about tones, 
fundamentals, and overtones is the 
assumption that a repetitive sound can 
be decomposed into a series of sine 
waves of different frequencies. That 
assumption, originally due to Daniel 
Bernoulli, is the basis of all Fourier 
analysis: to approximate an arbitrary, 
periodic function with a series of sine 
(or cosine) waves. Why would we 
want to do that? Mainly, because we 
know a lot about how to manipulate 
sine waves — we can easily take their 
derivatives and integrals, apply trig 
identities to them, and so on. We prob- 
ably can't do that with functions of 
arbitrary complexity. The Fourier 
series is yet another example of an 
approach I've stressed in these pages 
many times before: If you can't solve 
the problem you have, turn it into one 



that you can. 

Once we have a signal in its Fourier 
representation, we can also do tricks 
with it — ^the mathematical equivalent 
of tone controls. We can emphasize the 
high or low frequencies, filter out 
unwanted frequencies, and so on. Tone 
controls do that using analog circuitry, 
and we can do the same thing digitally, 
using a computer or digital signal 
processor. Done properly, tfie digital 
approach can be much more effective 
than any analog circuit, and is the 
foundation of digital signal processing. 

When J.B. Fourier first published 
his technique in 1822, the concept 
must have seem< d wonderful and mys- 
terious and perh ips even a bit risky to 
his colleagues. 1 . at thanks to that ubiq- 
uitous stereo system, to us it seems as 
natural as brea ling. You know the 
method works b cause you can hear it 
wifli your ears. 
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tions (/>n(x). The subscript n is there 
because we must have more than one 
such function, or else ^x) would have 
to be f(x). We can write the approxi- 
mation as a linear combination of the 
base fiinctions: 

fix) = oq^ q{x) + a\(t> |(x) + 

a2(t> 2{x) + --- + a„4)„{x) + --- (1) 

(The trailing ellipsis is meant to 
indicate that the series might be infi- 
nite.) Expressed as a summation. 
Equation 1 becomes: 



range. Much more likely, we must look 
at it only over some finite range, say r 
to s. Now let's do something complete- 
ly crazy. With no justification at all, 
we'll multiply f(x) by one of our func- 
tions, ^mW- Then we'll int^rate over 
tiie range to get: 



/ <t> m(.x)f(x)dx = 
r 

S CO 

/ ^^^n<P mix)^ nix) dx 



aaiA s), ! le whole term becomes a sim- 
ple con tant, computable beforehand. 
Let's ci II it Pn„. Obviously, we have 
more tl in one such constant — in fact, 
a potei tial infinity of them, one for 
each cr mbination of m and n. But con- 
ceptua y, we can organize them into a 
matrix Also, we can define a vector u 
such tl at: 



n=0 



«„= <l>„ix)fix)tk 



(5) 



/(x) = ^a„</>„U) 



(2) 



n=0 



Of course, we've seen plenty of 
series before, usually in the guise of 
poww series. You can now see that the 
form in Equation 2 is much more gen- 
eral; a power series is merely a special 
case of this more general form, where: 



We can interchange the order of 
summation and integration, which 
amounts to integrating term by term: 



/ ^ mix)fix)dx - 

r 

j (p mix)4> „ix)dx 



The e sments of this vector, of course, 
do cc itain fl[x), so we must perform 
the in cgrals for each case, but fbsy are 
straig itforward enough. 

W> can iu>w write Equatitm 4 in its 
equi\ almt niatrbc form: 



u = I'a 



(6) 



n=0 



(4) 



(l>nix) = X" 



(3) 



Our next challenge is to discover the 
right values to plug in for those coeffi- 
cients a„. At first glance, Has may seem 
an impossible task, but as you will see, 
it can be done even with the most gen- 
eral assumptions about the functions 
f„(x). 

First, we must recognize that we are 
extremely unlikely to be able to 
approximate f(x) over an infinite 



Well! It certainly doesn't seem lliat 

we've made our lives any easier — ^this 
is a messy-looking equation, contain- 
ing not just one, but both of those 
dread symbols, the integration and 
summation symbols. But take another 
look. The integral term in brackets on 
the right doesn't contain f(x), only the 
functions 0„(x), which are known and 
unchanging. So once we define our 
function set (and the range limits, r 



Ignc e for the moment some embar- 
rass ng technical details, such as the 
fact that we chose our indexing in 
Equ ition 4 to run firom zero instead of 
one and that a, a, and P are all poten- 
tial ! V infinite in size. As my major pro- 
fes; or used to say, "conceptually 
there's no problem." 

( onceptually, then, we can invert 
Eq lation 6 to find the desired coefli- 
cients, as shown: 

a = P-»u (7) 
As you can see, we really didn't com- 
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When we evahniie 
the integrals at 
their limits, all the 
sine functions 
cancel and we're 
left with lero. 

plicate the problem by multiplying and 
integrating, we solved it Conceptually, 
that is. You might still be wondering, 
however, whatever possessed us to 
take such a seemingly bizarre step. The 
short but somewhat unsatisfying 
answer is, "it works." The best answer 
I can give you is that we must some- 
how get into the solution, the behavior 
of f(x), not only at specific points but 
over the entire range of interest. An 
integral is the only way to do that. As 
for multiplying by 0m(x), we must also 
somehow isolate the eflFect of one of 
the coefficients from the total. Because 
squaring a function is sure to enipha- 
size its effect, the product seems to 
draw the influence of that one function 
out of the noise. The trick will become 
much clearer in the next step. 

You see, so far we haven't said any- 
thing about the orthogonality of the 
functions; in &ct, we haven't evai 
defined the term. But look again at Ike 
form of the matrix P. Since: 



Pmn=j4> mW0 nix)dx 



(8) 



the matrix is symmetric — = p^. 
Setting m = n gives us the diagonal 



(9) 



Because the function is squared, these 
terms must all be positive, and can be 



range (a pretty unlikely event). But 
what about all those off-diagonal 
terms? Wouldn't it be nice if, by some 
miracle, they all turned out to be zero? 
In that case, the matrix would be diag- 
onal, and its inversion would be trivial. 
In the sununation of Equation 4, only 
one term would be non-zero, so the 
eipaMon would become, simply: 



/ ^ „ix)fix)dx = a„p„, 



or. 



On = 



(10) 



where u^ is defined as in Equation 5. 
Our seemingly impossible task has 
suddenly gotten very simple indeed. 
Unless we're awfully silly in choosing 
the functions ^„(x), the integrals defin- 
ing will be available analytically, 
so we need only compute the values of 
u„ to get our desired coefficients. 

A set of functions ^„(x) having the 
property that the integrals of their 
products is zero, for all m and n except 
m=n, is called an orthogonal set. There 
still remains the question of whether 
we're asking the impossible in seeking 
such a set. After all, we're talking an 
infinite number of such functions; 
what are the chances that the integral 
be zero for all those infinite products? 
Fortunately, we can easily show that at 
least one set of such functions exists: 
the sine and cosine functions. Coiisider 
the integral: 



/ = jsinmxsmnxdx 

—n 

A trig identity gives: 



(11) 



The integrand is now simply the differ- 
ence of two cosine functions, and the 
int^ral is straightforward: 

' = ^5i^[sin(n-m)xl!„- 



2(n+m) 



sin(« + m)xY_ 



Now, the sine function is an odd firnc- 
tion. Not odd as in peculiar, but odd in 
the sense that: 

sin (-x) = -si/1 X 

When w e evaluate the integrals at 
their limits, all the sine functions can- 
cel and we're lefl with zero, as we'd 
hoped, foi- all cases except the <me 
where m = n. That case must be treated 
separately: 



/ = j sin- tixdx - Jid - co%2nx)dbi 
-It -n 

=k:=¥=- (13) 



We get similar results for the case of 

cosine functions and a combination of 
one sine and one cosine function. We 
have one more special case to consid- 
er, and that's the one for cos nx, where 
n = 0. This, of course, yields the con- 
stant "function" 1. The integral of (he 
product of this function with the sine 
function is, of course: 



/ = J sinnxdx 



which is clearly zero and the same for 
the cosine function. There remains 
only the integral of the function multi- 
plied by itself, which is, of coiose: 



sinasin/j = i[cos(a -b)- cos(a + b)] / = j ^ = n - {-n) = In 



(14) 



Thus: 

n 

I=j j [cos(n - m)x - cos(n + m)x]dx 

(12) 



So \' liat does all this tell us? Simply 
thai ,t's really possible to define a set 
of 01 ihogonal functions because the set 
consisting of 1, sin nx, and cos nx is 
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just such a set, on the interval -ix..n. 
This set of functions is, of course, the 
basis for the Fourier series. Equations 
13 and 14 also give us the values of Pn„ 
for this function set 

A DIGRESSION 

We've seen that sine and 
cosine waves form an 
orthogonal set. Clearly, the 
terms in a power series are not or&og- 
onal functions, since: 

will not be zero over any interval for 
all m and n. However, take a look at 
the first two functions in such a series: 

W = 1 

r,(jt) = x 

These two functions are, in fact, 
orthogonal. If we evaluate the cross 
integral: 

1 

/ = /ro(x)r,(jc)dic 

-1 

we get: 

/ = x|!l = l-l = 

At least we have two polynomial func- 
tions that are mutually orthogonal over 
a specific interval. Can we find a third? 
We know that the term Xj by itself 

won't work, but what about a different 
quadratic in x? Yes, the function: 

r2(jt) = 2*2-1 (15) 

is orthogonal to both other functions 
over the specified interval. Sharp-eyed 
readers will spot the three Amctions as 
the first three Chebyshev polynomials, 
of which there is an infinite number. 
These polynomials do indeed compose 
another infinite set of orthogonal func- 
tions. Not surprisingly, many other sets 
have also heea defined, each best suit- 
ed for deali^ irittt a pittotaf Mad of 
problem. 



BACK TO FOURIER 

I've already told the punch line of 
my story by giving you the orthog- 
onal functions used in the Fourier 
series — ^they are the sine and cosine 
functions. Specifically, if a given func- 
tion is periodic with period 2;r, the 
Fourier series c(»Te^ndiiig to tiw 
fimction is: 

fix) = ^ + a\ cosjt + 02 cos2jc + 

+03 cos3x + bi sinx + 
+^2 sin 2x + &| sin 3jc + • • • 

at, more concisely: 
fix) = ^ - f ^ ^ {a„ cos nx + b„ sin nx) 

n-l 

The first, constant coefficient is usual- 
ly separated out and given the factor of 
one-half shown, so we may use a sin- 
gle definition for all the coefficients: 

a« = i j f(x)cosnxdx (n = 0,1,2...) 

(17) 

and: 

6- = i / fix)^mtadx (n = 1,2,3...) 

(18) 

The multipliers of these integrals come 

from the definition of p„„ in Equation 
9, and our integrals obtained in 
Equations 13 and 14. Hie extra factor 
of 1/2 in Equation 14 is the reason for 
the factor of two in the ag term. We 
have no bo term because, of course, 
sin(O) = 0, so the term is zero regard- 
less of the value of the coefficient. 

One of the great advantages of the 
Fourier series is that it can handle fiinc- 
tions that are poorly behaved. 
Specifically, it's valid for functions 
that have discontinuities in them. This 
is very important to mathematicians 
because most math techniques such as 
differentiation require functions to be 
continuous. Using the Fourier series. 



we replace a fimction with discontinu- 
ities by a summation of an infinite 
number of fimctions, each of which is 
itself continuous. That's a neat mathe- 
matical hat trick. 

When Fourier first introduced his 
series, there was much concern over 
yAnsa and if we could be sure the series 
converged. Fortunately, it was later 
shown by Dirichlet that it does con- 
verge as long as certain rather loose 
conditions are met. These are the 
Dirichlet conditions: 

• The iunction must have at most a 
finite number of discontintiities. 

• It must have a finite inind)^: of max- 
ima and minima. 

• These maxima and minima must be 
bounded. 

These conditions are easily met for all 
but pathological fiinctions, so for all 
practical applications, we can be 
assured that the series will indeed con- 
verge. Also, any integrals or deriva- 
tives of the series will converge to the 
integrals or derivatives of f(x). At dis- 
continuities, the series converges to the 
midpoint of the extreme values. These 
characteristics delighted the matfie- 
maticians because they let them tame 
otherwise unpleasant functions. 

SYMMETRIES 

We can sometimes make use 
of symmetry properties to 
simplify the Fourier series 
for a function. We can split the inte- 
grals of Equation 17 and 18 into two 
parts like this: 


ci„ I f(.x) cos nxdK + 

— jr 

jr 

j fix)cosnxdx ^ 


71 

= J fi-x)cosnxdx-i- 



n 

^ J fix)cosnxdx 
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or: 



<^ = [/W + fi-x)\c(Knxdx 



Similarly: 

bn = nj if^x) - fi-x)]s:innxdx 



(19) 



FIGURE 1 

Even and odd. 



(20) 



An even function is one that is sym- 
metric about tbe X = point, as shown 
in FiguDe la. For mdli a fimction: 




An even function 



f(-xhf(x) 

and the coefficients become: 

n 

a„=if f(.x)cosnxdx 



(21) 



(ev^ functions only) (22) 



&,=0 



Similarly, a function is odd or anti- 
symmetric if its left side is a negative 

version of the right one, as shown in 
Figure lb. For this case: 



f(-x)=-m 



(23) 



For odd functions, the equations for the 
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fl[x) = -l 



(-n^<x<0) (25) 



(odd functions only) (24) 



*» = t/ fix)smnxdx 



Either way, half the coefficients van- 
ish. This simplification works in our 
favor in two ways. Not only do we 
reduce the number of coefficients to 
compute, but also, when we evaluate 
the series, we can skip the computation 
of either all the sines or all the cosines. 
Needless to say, it pays to look for 
such symmetries. Bear in mind, also, 
that functions may be symmetric about 
points other than x = 0. We can still get 
the advantages of symmetries, diougji, 
by a change of variables. 

APPROXIMATIONS 

Theoretically, the Fourier series 
for a function doesn't just 
approximate that function, it is 
the function, exactly. In practice, we 
never really get a perfect fit to the 
function because we can't compute an 
infmite number of coefficients or add 
up an infinite number of sine waves of 
increasingly higher frequency. 
Fortunately, we don't have to. Despite 
the aspirations of the more rabid audio- 
philes, an amplifier doesn't really have 
to pass frequencies somewhere above 
gamma rays to sound good. A fi'equen- 
cy response beyond the limit of audi- 
bility is sufficient. In the same way, in 
many applications of Fourier series we 
can get by with approximating a func- 
tion by its Fourier series eiqtansion, 
truncated to some finite number of 
terms. 

At this point, most textbooks offer 
some examples to show how well (or 

poorly) a Fourier series fits a given 
function. I won't break with tradition; 
let's compute the series for a few sam- 
ple inputs. If we choose simple func- 
tions, we can get the integrals in closed 
form. Consider, first, the ubiquitous 
square wave, shown as an odd fiinction 
in Figure 2. For this fiinction: 



bn = ij sin/ijcdx 
Since the function is odd. Equation 24 o 
applies. The integrals become: Int^rating gives us: 
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f(x)=l 



i0^x<n) 



CIRCLE # 48 ON READER SERVICE CARD 



= ^[-cosior + cosO] 

which gives: 

*» = 7^(l-cosfwr) 



FIGURE 3 

19th order. 



(26) 



Now, cos mr has a value of 1 , if n is an 
even number, md -1, if it's odd. So we 
can write: 



b„=0 



(n odd) 

(n even) 



(27) 



The coefficients for the first few 
terms are shown in Table 1 . The graphs 
of the first three approximations 
(orders 1, 3, and 5) are also shown in 
Figure 2. As you can see, the first 
approximation is simply a sine wave 
with amplitude 1.27. As we might 
expect, the approximation gets better 
for increasing order. The 19th-order 
approximation, shown in Figure 3, is 
beginning to look quite good. This fig- 
ure gives us a good feel for how the 
Fourier series handles discontinuities. 
Sudden, discontinuous jimips in ampli- 
tude are replaced by steep but continu- 
ous changes, with a little overshoot and 
ringing nearby. As promised, the value 
at the discontinuity is the average of 
the two extreme values. The little rip- 
ples in the tops of the wave are also 
typical of Fourier approximations. 
Looking at the figure, it's easy to see 
that the effect of the next higher term 
would be to sraoo& out diose ripples 

TABLE 1 

Coefficients for square wave 



n 


b„ 


1 


1.2732395 


3 


0.4244132 


5 


0.2546479 


7 


0.1818914 


9 


0.1414711 


11 


0.115749 


13 


0.0979415 


15 


0.0848826 


17 


0.0748964 


19 , 


0.0670126 




(and add new, higher-fi'equency ones). 

From the graphs, we can make two 
important observations. First, 
Dirichlet's analysis guarantees that the 
series always mm^:^, but it doesn't 



promise how fast. In fact, the conver- 
gence of the Fourier series is typically 
slow, which means we will need many 
terms in the series to get a good 
approximation. Stated in audio terms. 
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we need to reproduce many harmonics 
to accurately represent the sound of an 
instrument, even a low-frequency one. 
We could have predicted this result 
fix)m Equation 27; the amplitude 
decreases only as 1/n. To get an accu- 
racy of one part in 1,000, we'd need 
500 trains. Obviously, we're not going 
to get the kind of part-per-million 
accuracy we're used to in other numer- 
ical approximations. Fortimately, we 
rarely need that kind of accuracy, and 
an approximation as good as the 19th- 
order fit of Figure 3 is often good 
enough. 

Second, if you compare the over- 
shoot in Figure 3 to diose in Figure 2, 
you get the impression that it's not 
really decreasing much — the over- 
shoots for tiiird, fifth, and 19tfa order 
seem almost the same. You might well 
ask, how many terms do I need before 
the overshoot starts decreasing? Itie 
answer it never willi 



Surprisingly, the overshoot has a 
constant amplitude for all but the low- 
est-order fits. However, for an approx- 
imation of any order, we can find an 
equation for die maxima and minima 
with a little calculus. The overshoot is 
the value of f(x) at the first maximum. 
If we evaluate this for various orders, 
we find that the height of this maxi- 
mum does not converge to zero, but to 
a value indq)endent of n, and roughly 
11% larger than the value of the origi- 
nal function f(x). This effect is called 
die Gibbs phenomenon. 

The phenomenon may seem to con- 
tradict the statement that the series 
converges to the original ftmction. 
Remember, this effect only occurs in 
the vicinity of a discontinuity, where 
all bets ms off anyway. At any point 
away from the dkcoBtinuity, even by 
an infmitesimal distance, the series 
behaves iimnally. In practice, the 
GiUis ^i^aamenon is a minor annoy- 



ance, nothing more, albeit an interest- 
ing and sui prising one. 

Let's try another function, the saw- 
tooth, which is shown in Figure 4. For 
this function: 

f(x) = x (28) 

This is another odd function, so flie 
coefficients are given by: 

an=0 

(29) 

b„=^ j xsinnxdx 



A table of integrals gives: 
j xsinxdx = sinj: - jccosjc 

This equation has the same form as our 
equation for b^ that has the factor, n, in 
Hhs sine functi<m. If we let y = mt, 
Equation 29 becoi^: 

^n=-^ Jysinydy 



Our formula then gives tis: 

^^Isiny-ycosylo" 

= ^[(sin»«r - RTrcosinr) - 

(sinO-OcosO)]' 
= ^(sinior - lurcosiur) 

Now, sin nTT = for all n, and cos wr is 
either -1 or 1. We get, finally: 

bn=ii-ir*' (30) 

Interestingly, the amplitude of th> 
coefficients goes again as 1/n, just a 
for the square wave. The difference 
that this series contains all harmonic 
not just the odd ones, and they alte 
nate in sign. The first few odd orders 
the approximation are shown in Figi 
4 and the 19th order in Figure 5. Aga 
note the overshoot (Gibbs phenor 
non) and ripple that are cbatacteri: 
of the approximation. 




PnMfilCE 

T Ethernet Support 
▼8051 Debugging 
▼XRAY+/pSOS Support 

cm TODAY! 1 -800-PROMICE 
(1-800-776-6423) 



Grammar Engine Inc. 

921 Eastwind Dr., Suite 122 
Westerville, OH 43081 
Phone 614/899-7878 
Fax 614/899-7888 




CHANGE OF VARIABLES 

By now, you're probably con- 
vinced that the Fourier series 
really works and provides a 
good approximation for any periodic 
function. So far, though, we've only 
shown functions that are periodic with 
period In and centered at zero. Real 
functions aren't always so accommo- 
dating. In general, a periodic fiinction 
of time has period T and, equivalently, 
a ftequea^ given by: 



FIGURE 4 

Sawtooth wave. 



f = 1/T 



(31) 



(Sorry about the duplicate use of the 
symbol f, but its use for both the fre- 
quency, f, and the function, f(t), are 
rafter traditional. When 1 refer to the 
function, I will always do so wifli the 
parentheses in place, so there should be 
no confusion.) 

We'll usually find it convenient to 
specify the function over the range 
0..T, rather than centered at zero. So 
we need a relationship between x and t 
such that: 

when X = -jT, t = 
when X = ff, t = T 



The transformation: 



(32) 



does the job. Revotmg this relation- 
ship gives: 



(33) 



Many texts show the Fwirier series and 
its formulas for the coefficients in 
terms of t rather than x. Personally, I 
find it easier to change the variable as 
needed externally and keep the formu- 
las in their simpler and more pristine 
state. 

IMPLEMENTATION 

So much for theory. Now, let's 
put the theory to work with some 
software. We'll need software of 
two kinds: one to compute the Fourier 
series, given x and the coefficients, and 




the other to compute the coefficients 
for a given function. The first is by far 
the easier. We'll assume that the eoef- 
ficients are given to us as vectors a and 
b (for once, ^ C ©Hvemtion of index- 
ing ftom zero woite l» w ^vant^e). 



A direct translation of Equation 16 
gives the code shown in Listing 1 . Sure 
enough, this code gives, the desired 
result and in only four lines of exe- 
cutable code. 
Unfortunately, this is ^so ttie worst 
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USTINGI 

Fourier first cut. 



FIGURE 5 

19th order. 



double fourier(double x, Vector la, 
Vector lb, int n){ 
doubiLe result - a[0]/2.0; 
for(iJit i=l; i<=n; i++) 
result += a[i]*cos(i*x) + 
b[i]*siJi(i*x); 
return result; 

way to compute tiie series. The prob- 
lem is that we end up calling the sine 
and cosine library functions 19 times 
each. But, you ask, how can I avoid 
that? The series depends on them. Ah, 
but they're not independent. Each of 
the angles is an integral multiple of the 
original angle, x. We can use this fact 
and some trig identities to reduce the 
number of function calls to two. 

A pair of trig i^itities ^ye& us the 
following: 




sin(a+^steae0S&+cosaan6 (34) 
cos(a +6)= cos a cos 6- sm a sin j> 

IM a » X, and b = nx. We can then set 

up a recursive method to conqnite the 

next pair of functions: 



8in((iri- 1 )x)=sinxcosnx +co8x»nnc 
co8((r+ 1 )x)>=cos xcosnx -sinordn/K 

(35) 

Since we only visit the functions for 
nx one time each, we can update sin nx 
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and cos nx in p^ice, like tiiis: 

tegip > sin_x * cos.nx cos_x * sln.nx; 
cos_nx = cos.x * cas.nx - slnjt • sinjn; 
sin_nx - temp; 

We should start with nx = 0, which 
gives sin_nx = 0, cos_nx = 1. Listing 2 
gives the improved code. It's not four 
lines long anymore, but it's still quite 
conqiact. More to the point, it's about 
as efficient as you can get. I suppose I 
should point out that we never use the 
value of b[0]. But it is easier and less 
confusing to carry that coefficient 
around, keeping the indexing as is and 
wasting one storage location, rather 
than trying to save that single word and 
obfoscating ttie code unnecessarily. 

COMPUmiG THE COERiCffillS 

In our two examples, the square wave 
and sawtooth, we were able to get ana- 
lytical expressions for tiie coefficients 



a„ and bn- That's because I chose ftinc- 
tions such that the integrals in 
Equations 17 and 18 were integrable in 
closed form. Clearly, we won't be able 
to do this in general, and a computer 
implementation must deal with the 
issue of how to compute the coeffi- 
cients. This is where the rubber meets 
the road because numerical integra- 
tions are slow, and we must integrate 
over the entire range -Jt..7t, once for 
ev«y coefficient. 

The first decision we must make is 
what kind of integration algorithm we 
should use. In my previous articles on 
integration algorithms, I covered all 
kinds of fancy schemes, including 
high-cwder difference methods like 
Adams-Moulton and single-step meth- 
ods like fourth-order Runge-Kutta. 
Wtcm tfie conclusions in those articles, 
we're tempted to say that the higher the 
order of integration we use, the better, 
beeaose we can use a larger step size. 



USTINGZ 

Fourier improved. 

dMble fourierCdoufale x, Vector ta, 
Vector ftb, int n){ 
i (toulxle siji_x ' sin(x); 
double cos.x = cos(x); 
double sinjjt = 0.0; 
double cos jx,^;^, , 
double result » a[()]/2.0; 
double teap; 

for(int i<=n; i++){ 
temp = sin.x ♦ cos^ix + cos.x 

* siji.ix; 
cosjjt = cos.x * cosJLx - siii.x * 
« sinjji; 

sinjx = teup; 
i,^ result ali]»cosJji + 

b[i]*sifljt-'v^ <j- ;| 

return result; 

■ } WW 

We'd be wrong, and therein lies a tale. 
Many years ago, I was trying to inte- 
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grate a complicated function of an 
angle, 0, over the entire range, 0..2n. I 
was looking for an analytical solution, 
but being unable to find one, I resorted 
to numerical integration. For simplici- 
ty, I chose a trapezoidal integration 
(second order) with lots of steps, I 
think about 1,000. 1 got an answer, but 
wasn't sure how accurate it was, so I 
performed the usual sanity check: 
increase the step size (reduce the num- 
ber of steps) and try again. I got the 
same answer. Not just nearly the same; 
exactly the same. "Hm," says I, "I 
guess I didn't need as much accuracy 
as I thought." So 1 doubled the step 
size again, and got the same resuh. 

Aftet I got down to only 40 or 50 
st^, I was beginning to worry. The 
result should have started showing the 
effects of roundoff error long ago, but 
it was still hanging in there, rock solid. 
I felt that something must be wrong 
Willi my algorithm; it was going to give 



the same result regardless of the step 
size. Eventually the roundoff error did 
start to appear, but only after the num- 
ber of steps was only 10 or so. This 
result seemed to defy intuition to such 
an extent that it was quite awhile 
before I was able to accept the answer 
as valid. Finally, however, I figured 
out what was going on. 

If we're going to integrate a function 
analytically, we must furst have a table 
of its values or at least be able to gen- 
erate such a table. The integration pro- 
ceeds in discrete steps with some step 
size, h. All numerical integration algo- 
rithms are based on the notion of fitting 
a polynomial of some degree through 
these tabular points. Integrating this 
polynomial for one step gives us an 
integral that's a linear combination of 
adjacent tabular values: 

(36) 



where f.„ stands for the tabular value 
of f(x) n steps behind the current one. 
The values of the coefficients c„ 
depmi, of course, on the order of the 
intention scheme. But we can be sure 
of one thing: They must add up to 
unity, because if the function is con- 
stant, the integral must be equivalent to 
hfo- 

Now let's apply Equation 36 over 
the entire range of the function. We'll 
end up summing all the individual dis- 
creet integrals to get: 

/ = ^y^(co/o + ci/_i + C2/-2 + 

•••c„/_„) 

on 

(37) 
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In a general integration of a function 
that goes on forever, the sums in 
Equation 37 will all be different 
because they cover different ranges of 
the tabular points. But if we're dealing 
wMi a p^oiHc fiuietion, tiie values of 
fj start repeating each time we cross the 
limiting value (In, in our case). All the 
sums are the same, representing a sim- 
ple sum of all entries in the table. We 
can let: 

S = Y^fi (38) 

where we understand that this is to 
include all the tabular points once and 
only once. Then Equation 37 becomes, 

1 =5 h(€Q$ + cjS + cjS + • • • + c„iS) 
»ii|ei) + ci + C2 * — + c„) 

But we've already said that the coeffi- 
cients must add up to Mnity, so we 
finally get: 

l = hS = hJ2fi (39) 
j 

This is a vitally important result. The 
integration method is seemingly first 
order, corresponding to a simple sum 
of all the tabular points. But what 
we've seen is that it's equivalent to fit- 
ting a polynomial of second order, or 
third order, or, in fact, any order up 
through n, the number of points, 
through the tabular data. In other 
words, the integration method may as 
well be nth order, which would give us 
an incredibly accurate algorithm. 
That's why my simple-minded 
approach worked so well. 

We can use this same approach to 
get the integrals we need for our 
Fourier coefficients. No fancy integra- 
tion schemes are needed; m fact, a 
fancy scheme would be wasted. One 
word of caution: If you're tabulating a 
function of an angle that runs over a 
range something like 0..27r, each tabu- 
lar point must be included once and 
only once. It's tempting to include both 
the value at zero and the value at In to 
make the function look pretty in a doc- 
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ument. But to do so would result in that 
point being added twice. Avoid the 
temptation. 

The simple addition of the data 
points works fine when a function is 
continuous, but what do we do about 
discontinuities? The temptation is to 
include two points at the same value of 
X, but that would mess up our counting 
scheme and integration formulas. The 
key to the right approach is to remem- 
ber how the series behaves at a discon- 
tinuity: It gives a value equal to the 
mean of ffie Wo extreme values. If we 
provide the same value at a discontinu- 
ity, we'll get the right answer. 

Of coiffse, ^t approach oniy works 
if the discontinuity lies on a mesh point 
of the table. If it doesn't, we can only 
^ve tiie vsdiies fiar the two bracketmg 
points and hope for the best. 
Incidentally, my statements regarding 
the hig^ acG^na^, even when we have 
few data points, don't apply for this 
case. It should be clear that fitting a 
polynomial through discontinuous 
points is tricky, at best. When a fiinc- 
tion is discontinuous, there's no substi- 
tute for a lot of data points. 

THE SUM OF IT ALL 

Let's see what our integrals look 
like when replaced by summa- 
tions. Equation 17 gives (with a 
change ( 



fly = i / f(.x) COS jxdx 



As the coefficients (and the Fourier 

series) are defined, the index is sup- 
posed to run to inllnity, but uf counsc 
we must limit it to some finite number, 
say m. 

Now, suppose that f(x) is tabulated 
with a total of n points, and let i be the 

index into the table. The st^ size or 
granularity of the table is: 

h = ^ (40) 

The first point, for i = 0, corresponds to 
the point x = -tt. The second, for i = 1 , 
corresponds to x = -;r + h, and so on. In 



general: 
fi^f{-n + ih) 

Our integral is equi^ent to: 

aj = cos;(-7r -I- iTi) 

j 

Substituting for h from Equation 40 
gives us the multiplying coefficient in 
tenns of n: 



The formula for bj is similar: 
bj=iJ2^iSmji-n+ih) (42) 

A straightforward implementation 
would sum over all i, as in die follow- 
ing pseudocode: 

for j = to n-l Idop 
alj] = bij] = 
for i = to n-l loop 
X = -pi + i*h 
a[j] += f[i] * cos(j * x) 
b[j] += f [i] ♦ sin(j » x) 
end loop 
end loop 
a = 2*a/iti 
b » 2*b/B 

The problem with this approach is 
that we must have access to the entire 
array f [] each time through the inner 
loop. That may not be a problem if the 
array has reasonable size, but if it is so 
large that it must be read from a file, 
we would have to read that file m 
times. We can improve flie algorithm 
greatly by exchanging the order of 
summation: 

a » b « 

for i = to n-l loop '' 
X = -pi + i*h 
for j = to n-l loop 
a[j] += f[i] * cos(j * x) 
bIj] += f [i] * sin(j * x) 
end loop 
end loop 
a = 2*a/n 
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b»2*b/i 

Now each table entry is visited only 
once in sequential order. 

Of course, we have the same effi- 
ciency problem as before. As the algo- 
ridun is written, we end up calling the 
sine and cosine fianctions mn times 
each — a prohibitively costly approach. 
But we can use Hie same trick as 
before, our double-angle formulas, to 
compute each trig function from the 
previous one and reduce the number of 
functions calls to two. This time, we 
can apply the trick twice, once for each 
loop. 

Sharp-eyed readers will recognize 
this process as strength reduction, a 
technique for loop optimization that 
we studied several months ago. In this 
technique, we replace a multiplication 
involving a loop index by an equiva- 
lent addition. Huis, for exao^Ie, the 
code: 

for i = to n-1 loop 

X ■ -pi + i«h 
end loop 

becomes: 

X = -pi 

for i = to n-1 loop 

X *= h 
sid loop 

The only difference is that now we are 
not only applying strength reduction to 
the multiplication operator, but to trig 
fiinctions as well, resulting in a huge 
reduction in overhead. 

IS IT rn YET? 

Practical code for the algorithm is 
shown in Listing 3. The calling 
convention is simply: 

fMiifitlf, a, t>); 

where the array f contains the input 
data, and a and b are the ou^ut coeffi- 
cients. 

Both f and the coefficient vectors 
aie declffl«d as class Vector. As I've 



defined this class, each vector carries a 
length parameter, so no additional data 
in the way of size parameters need be 
supplied. The order of the fit is deter- 
mined by the ranks of vectors a and b, 
which should be of equal size. If they 
aren't, the algorithm will use the rank 
of the shorter vector. Remember that 
the order of the fit includes the con- 
stant tettn tM, so the size of the coef- 
ficient vectors should be the desired 
order of the fit, plus one. These vectors 
must be deeiaiei itad sized in the call- 
ing program. 

As you em see, I've implemented 
fourfitC) M lUvmi^ the tabular data 
defining the function is stored in an 
array. However, since each data point 
is needed only once, you can easily 
modify the code so each element of <(> 
is read from a data file instead. As a 
matter of ^kA, W f is treated as a virtu- 
al memory array, and if such an array 
is defined as a C++ class, no code 
changes will be required at all, other 
than changing the type of <p. 

Wilk this algorithm, we are getting 
vetf dose to die famous Fast Fourier 
Transform. But no, we're not quite 
there, yet. In the first place, remember 
that me'te dealing witii the Fourier 
series here, good only for periodic 
fimctions, and not the Fourier trans- 
form. The latter is a integral relation- 
ship. But since we've already seen how 
integrals get implemented as summa- 
tions anyway, we're cIosct than you 
might think. Because our strength- 
reduction technique is quite similar to 
that used in the FFT, we get very good 
efficiency. In my tests on a 486-based 
PC, the algorithm processes a 19th- 
order, 1,024-point transform in well 
under a second. I call it the "half-fast 
Fourier transform (FFT/2)." 

In one sense, the algorithm is actual- 
ly more general than the real FFT. 
Hi-iMu.s'o of ihc way the loops are 
implL-inciiici.1, iluii algorillim ivijuircs 
that n, tlie number of data points, be a 
power of two. The FFT/2 algoriflim 
does not — any number will do, as long 
as it is greater than the number of 
coefUcients. 



usimes 

Fourier coefficients 



•oid fourfit(Vector U, Vector la, 
Vector U)){ 
-int n * f.rankOr^ 
lilt m = niin(a.rank(), b.rankO); 
doulile fcj, fsj, tenp; 
doubUe coeff » 2.0/n; 
double h = coeff *pi; 
double s = sin(h); 
^double c = cosdiJi;" ' 
double ci = -1.0; 
double si = 0.0; 
'^t^:-* 0.0; 

•> = 0-0; 

for(iJit i=0; i < n; i++){ 

fcj » f[i] ♦ ci; 
fsj » f m » si; 
for(ijit j"l; >j .< n; }**){ 

a[J] += fcj; 

bU]*=fsj; 

tenp » fsj » ci + f cj * si; 
fcj = f cj * ci - fsj * si; 
fsj = temp; 

} 

tenp » si * c + ci * s; 
ci » ci * c - si » s; 
si s tenp; 

} 

a *= coeff; 
b *= coeff; 

We still have one more step to go, 
extending our result using the Fourier 
series (a summation) to the limiting 
case of an integral. We'll do that next 
time, and that will give us a full-up 
FFT implementation. We'll also look 
at a complex-variable version of all the 
functions. Stay tuned. 

./cicli W. Crenshaw is a staff scientist at 
In Vivo Research in Orlando, ft. Ik 
diii much early work in the space pro- 
iiriini and has developed numcr.nis 
analysis and real-time programs and 
holds a PhD in physics from Auburn 

University. Crenshaw enjoys contact 
and can be reached via e-mail at 
72325. 1327@compuserve.com. 
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Our prodigal mathemagician returns to the 

realm of Fourier analysis, transforms, and 
other matters of infrequent interest. 



In the first part of this series, we 
considered the whole concept 
of approximating a function by 
a collection of o&er functions.' 
In other words, given a function 
f(x) defined over some interval, we'd 
like to i^lace it by a smes: 

Xx) = oo^ o(^) + fli^ 
+a2(t> 2(x) + — + a„(l)„{x) + -- (1) 

You're all fij^Uisr with &e notion of a 
power series — a special case where the 
functions 0(x) are simple powers of x. 

Why would we want to replace a 
single function with a whole (poten- 
tially infinite) set of them? The usual 
reason is that we know something 
about the new fiinctions, such as how 
to integrate and differentiate them, that 
we didn't know about the old one. This 
is yet another case of the time-honored 
principle you've seen before in these 
pages: If you can't solve a problem as 
it stands, turn it into oae you know 
how to solve. 

To make the transformation work, 
you must know how to compute the 
coefficients a„. We looked at a few 
ways to do that and discovered that the 
process is particularly easy if &e fimc- 
tions are orthogonal, which means &at 
they obe^ the relatimisfaip: 



/ 



(p „(x)<t> „(x)dx - for all m,n 

such that m ^ n (2) 



which the original function was 

defmed. 

For convenience, we often scale the 
functions so the integral is unity vAien 
m = n, but that isn't really necessary. 

It might seem at first glance that 
Equation 2 gives us an awfully strin- 
gent requirement, and you'd be right to 
wonder if finding such a set of fimc- 
tions is even possible. It is. In feet, a 
set of sine and cosine functions is such 
a set (though far firom the only one), if 
we choose the range and firequendes 
propffl-ly. Specifically: 



j sia mx sin nxdx 
-It 



(m ^ n) 
{m = n) 



im^n) 
2ji (m=«=0) 



j cosmxcosnxdx = 

-X 

/ smmxcosnxdx = \ 



) 



where r and s represent the range over 



(3) 

Notice that when m and n are both 
zero, it doesn't matter what x is; the 
cosine is constant (equalling one), so 
the integral is a littie different. The 
extra two in the result causes us some 
bother for this case, but we live with it. 

The notion of using sine and cosine 
functions to approximate a fvmction is 
attributed to Jean-Baptiste Fourier, so 
naturally the subject is called Fourier 
analysis. Last time, we saw that the 
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Fourier series is de^aed to be: 
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fix) = ■^+^^('^n cosiw + b„ sianx) 
«=i (4) 



(There's that pesky two.) The coefB- 
cieots and are defined by: 

On = i / /WcosiBBfa (n = 0,1,2...) 

(5) 

and: 

b„ = ^ j f{x)sinnxdx (n = 1,2,3...) 

(6) 

The idea of evaluating a potentially 
very large number of integrals to get 
the coefficients does not appeal and 
suggests very slow computations. But 
we can speed up the process markedly 
by recognizing that the integrals can be 
replaced by simple summations, and I 
showed you the justification for that. 
Finally, we learned that we can avoid 
tons of caUs to the sine and cosine 
functions by xmng ibe doubie-ai^e 
fcnmulas: 



sin(a + b)=siaacosb+aKaanb 
cos(a + cos a cos 6 -sin a sin 6 



(7) 



This trick lets us reduce a large nimi- 
ber of fimction calls to a mere two. 
Drawing on that technique, I gave yoU: 
the code for two functions that com- 
pose what I call, with tongue firmly in 
cheek, the "half-fast" Fourier trans- 
form (HFFT). The code is reproduced 
in Listing 1 . Function f ourier evaluates 
the Fourier series (Equation 4) from 
the coefficients, and fourfit computes 
the coefficients fhim an input function 
f (x), given in tabular form. 

Sharp-eyed readers will notice slight 
differences between this code and the 
equivalent listings fi-om the last install- 
ment. First, I've removed the number 
of data points fmm the argument list of 
fourier. This number can be inferred 
from the rank of the vector defining 



month, we'll 
look at ways to 
further Improve 
the efficiency of 
the MNpitation 
and use it in 
different 



fl[x). Second, I've replaced an error 
message testing the ranks of a and b 
with a C assertion. I've renamed sever- 
si vn^sMes in fourfit to make them 
more consistent and less Fortran-like. 
The changes don't modify the algo- 
rithms, just the readability. Last but fer 
fi-om least, I corrected a serious error in 
fourier — a sign error that made all the 
results have the wrong sign. 

This month, we'll look at ways to 
further improve the efficiency of the 
computation and use it in different 
ways. This will lead us to the discrete 
Fourier transform (DFT). Then, we'll 
look at what happens when we extend 
the concept of a series (a simi of terms) 
to an integral. This will give us the 
Fourier transform. Unlike the. Fourier 
series, the Fourier transform allows us 
to extend the concepts into the realm of 
nonperiodic fiuKtioss. 

SERIES OR IITTEGRAL? 

That last statement calls for some 
clarification. From some of the 
responses I've gotten to Part 1 
of this series, I've realized that many 
people are confused over the differ- 
ence between the DFT and the FFT and 
between the Fourier series and the 
Fourier integral. Here's the difference 
in a nutshell: The Fourier series works 
only for periodic functions. Because all 
of its terms contain factors like sin nx. 
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S'"'''tefflp = sin.x * cos.h + cos.x * 

. COS.X = cos.x * cos.h - sin.x * 
sin.h; 




the series is going to repeat itself every 
2jt radians, no matter what. Your input 
function f(x) may not be repetitive; it 
may be only defined on the interval 
7r..7r (or some ^propriately scaled 
equivalent) and be tt^sfiQed outside 
that range. 

Regardless of how the fimction f(x) 
behaves outside flris range, the approx- 
imation embodied in the Fourier series 
is goir^ to repeat itself at 2ji intervals. 
To the mathematics, fise vshtss -it and 
+n are the same point. The values for 
your fimction at these endpoints had 
better be equal as well. You can apply 
Fourier analysis to fimctions in which 
the values at -n and +n are different, 
and many people have. But you may 
not like the results because the mathe- 
matics tries to wrap the function 
aioxmd at those endpoints, and it will 
treat any differences between the two 
ends as a jump discontinuity. You saw 
in Part 1 tiiat ihe Foinier series handles 
discontinuities well. But you've also 
learned to expect "ringing" and edge 
transients smdk 8$ ClM>*s phenomenon 
in the vicinity of fliose discontinuities. 
If you tiy to use the Fourier series to 
approximate a segment of a general, 
nonperiodic function, you will see 
those effects, and the series will not 
give you a good approximation to the 
function, even within lie range over 
which it's defmed. 

The Fourier integral, on tilie olher 
hand, is another matter altogether. In 
general, it will contain components of 
all frequencies, and they won't all be 
commensurate. So, there will be no 
value of x for which the function 
repeats itself Fourier ttansfoims ^ply 
to functiom that are nonpeiodic, 



whereas the Fourier series caimot. 
That's the fundamental difference 
between the two. We'll go into all this 
in detail, but the distinction between 
the two cases will help clear some of 
the fog in what follows. For now. we'll 
continue with the case of the Fourier 
series, not the integral. 

Between the DFT and the FFT, there 
is no difference in fimctionality. If both 
are done correctly, the two algorithms 
should yield identical results. The dif- 
ference is internal and, as the name 
implies, comes fi-om seeking ways to 
speed up the DFT. Make no mistake, 
we're not talking about minor 
improvements in efficiency here. The 
FFT is not just faster, it's profoundly 
faster because as the number n of data 
points is increased. Has time required to 
perform the DFT increases by 
order(n2), while that of the FFT 
increases only by order(n log n). For 
large n, that difference can be stagger- 
ing. My cherished "half-fast Fourier 
transform" is really a DFT, so perhaps 
even that w^dm^cal name is fo' too 
optimistic. 

HOW MANY POINTS? 

Mathematically, the "input" 
fimction f(x) is a continuous 
function, defined over some 
range that we've taken here to be -n to 
+7r. Also mathematically, the ntmiber 
of coefficients is infinite. The series 
never quits except for very special 
functions. In the real world, of course, 
things look quite di£G»ent We can 
expect to be given the function not as 
an analytical one, but as a table of 
points. (If we had an analytical expres- 
sion for f(x), we'd use it). What's 
more, we can expect to get not only a 
table, but a table based on evenly- 
spaced divisions in x; the whole notion 
of using the double-an^e formulas 
depends on this. Finally, we must rec- 
ognize that we can't carry aroimd an 
infinite ntmiber of coefficients, so we 
end up truncating the series at some 
order. This step means that our series 
no longer precisely represents the input 
fimcti^ but is an ^roximatton to it 
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Let's call the number of input points 
n and the order of Ifae i^jproxiraation, 
or the highest frequency we retain, m. 
(Careful — these are not the same m 
and n as in Equations 1 tlm>u^ 6. The 
usage is conventional and the functions 
in Listing 1 are consistent with it.) 

When I first wrote &e functions 
fourier and fourfit, it was to use the 
concept of a series approximation pre- 
cisely as Fourier intended it to appton- 
imate a periodic ftmction f(x), given a 
single value of x. I was building 
iymonic simuli^cBis at die time and 
wanted to approximate things like 
aerodynamic effects. Since the angles 
associated witii aerodynamic effects 
are, after all, angles and inherently 
periodic, such effects are natural candi- 
dates for tills kind of sqiproximation. 
The idea was to execute function four- 
fit offline, using it only once to give 
Hbs nee^ coefficioits. These would 
then be appropriately truncated and 
stuffed . into the simulation as data. 
During a cranputation, function fouri- 
er would be called as needed for spe- 
cific values of x to give me the approx- 
anated forces and torques. Keep these 
functions in your toolkit because they 
are very usefiil for this kind of applica- 
tion, and you won't find ftmction 
fourier in any textbook I've ever seen. 
Its use to evaluate the series at a single 
point is rather special 

When using these routines for their 
intended purpose, I wanted as accurate 
a representation as possible witiiin per- 
formance limits. Since function fourfit 
would typically be run only once or a 
few times for each new simulation, I 
really didn't care how long it took to 
execute. If you'd asked me how many 
points I wanted in tiie teble for f(x), I 
would surely have smd, "As many as 
possiWe." 

In the execution of the simulation, 
on flie other hand, time was critical, so 
I wanted to use as few terms in the 
series as possible. Typically, I'd find 
that fifth or seventh order was ade- 
quate. If you'd asked how large m, the 
ordea:, should be, I'd say, "Something 
more than scvol" I needed to see only 



If we're doing 
Image 

enhancement or 
data cMpMekMi, 
we don't want to 
see oni point 
of tbe transformed 
image, but the 
wMe Hilng. 

enough terms to see how the coeffi- 
cients were decreasing with order, so 
I'd know where to sraies. 
In my case, n was much, mwl pei^ 
than m. 

You must understafii, 
this isn't the way most people use 
Fourier analysis. In my article, "The 
Alphabet fitm S to E,*^ we saw lhat, 
while transforms such as the Laplace 
transform and z-transform can be used 
in both diieetions, control S3rstems 
engineers often stop halfway. Except 
when Laplace transforms are being 
used to solve l^gH«atiaI equaticms, we 
find little need to reconstruct the fiinc- 
tion we started with. Control systems 
engineers an^ uraot to see the transfer 
function in s or z space because they 
can infer fi'om the nature of that fimc- 
tion what the behavior of their S3rstem 
is going to be. They don't need or want 
the inverse transform. 

It's the same way wifli Fourira- trans- 
forms. Commimications engineras and 
others who use such transforms often 
don't want to reconstruct the fimcticm 
f(x). After all, they already had it once. 
What they're interested in is the fre- 
quency and phase content of the signal. 
For that purpose, they need to 2see as 



much of the spectrum as possible. If 
ycm adc tiiose engineers how many 
co«sBBcients they need to see, they're 
likely to say, "All I can get." In this 
case, that term represents a specific 
number. According to Shannon's theo- 
rem, when a given signal is sampled at 
a cralain fiequraicy, no information 
contained in that signal can be 
retrieved if it's more than half the sam- 
pling fi«qu«icy. Producing a tabular 
representation of a continuous ftmction 
with, say, n points is equivalent to sam- 
pling at a firequency of n times ihe base 
frequency. Thus m can be no larger 
than °/2. Since we have both sine and 
cosine fimctions, we'll get out n coeffi- 
cients in total. And the folks who are 
looking at fi-equency spectra will want 
to see every one of them. 

A similar situation holds in the 
opposite direction. If we're doing 
someOung like im^ie enhancement or 
data compression, we don't want to see 
just one point of the transformed 
image, but the whole thing. These con- 
siderations represent the major differ- 
ence between functions fourier and 
fourfit as tiiey were intended to be 
used and the way most applications use 
the same formulas. In my case, the for- 
ward trmsform (fourfit) was a many- 
to-few transformation, while the 
reverse transform (the series evalua- 
tion, fourier) was a few-to-one. Most 
modem applications involve a many- 
to-many transformation. In fact, the 
number of output data itrans is almost 
surely the same as tbe number of inpat 
data items. 

The form of fourfit already fits this 
usage; we need only make the vectors 
the correct length (the coefficient vec- 
tors should be half as large as the ftmc- 
tion vector). We could get the other 
half (the vector version of fourier) 
simply by calling fourier n times, but 
that would hardly be efficient. But if 
you look at the form of fourier, you'll 
see that we don't actually use the value 
of X itself, but rather its trig functions. 
We can compute these fiinctions in 
recurrence relations, just as we have in 
fourfit. The end results looks awfiilly 
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familiar; it has almost exactly the same 
form as fourfit. The routine is shown 
in Listing 2. A side-by-side compari- 
son of inv.f our and fourfit shows that 
they're almost the same. This, of 
course, is as it should be; we saw in 
•The Alphabet from S to' Z" tiiat the 
inverse Laplace transform looks 
almost like the forward one. The rou- 
tines will look evea more alike when 
we take the next step, which is to con- 
vert them to use complex arithmetic. 

IHE COMPLEX SERIES 

As any electrical engineer 
knows, the relationship 
between trig fimctions and 
complex numbers is very close indeed. 
Specifically, a pair of identities gives: 



Equation 4 gives: 

fix) = ^ + 2^ + - 



ib. 



Collecting terms giv^ 

'a„ ib, 



Stnjc = f-(e" 



(8) 



■si 

We c^ iNdte efBite m 

00 00 

/W« Co + + ^ 4e-*« 



Substituting these relationships into ^ere: 
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(n>0) (11) 



(n<0) 



In the third sum of Equation 10, let's 
r^lace n by -&. Thesa Hie Stan 
becomes: 

B=-00 

Equation 10 now tsdces the foim: 

—1 00 

fix) = +52<i„«*«+Co+^c„e*« 

«=i (12) 



The first simi takes us over all the 
negative integers and the second one 
all the positive ones. The constant term 
takes care of n = 0. What's more, the 
exponential in each element is identi- 
cal. This holds for the constant term as 
well, since e'"" = 1 when n = 0. We can 
make the symmetry complete by let- 
tii^ d„ = c^,, implyii^ Hat: 



2 ^ 



ibn 
2 



(13) 



If a„ and b„ are real, as we've 
assumed, the for negative n are sim- 
ply the complex conjugates of their 
positive counterparts. 

Our Fourier series now talKS the 
very simple form: 



(14) 



In this representation, even that 
pesky power of two for the constant 

term has gone away. Equation 14 
has a certain elegance about it, don't 
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USTING 2 

Generating a result vector. 

//' Regenerate Fourier function 
, vcui inv::four,it^-i.i aa 

doublfi^osj);, sin. ix 
idouliLe teoip, residt: 
int'n;_^=?f;rank{); 
mt n = a;rank(): 
assertCa.rankO = o.rankO) 
double hr - Inijni. ' 
double cos.h = cos(h); 
aoueiB. sjji.n = san(h); 






for{jjit j=l; j<iii; j++){ 


1 


tes^ =-sirT X * cos IX * coijf. * ■ 






a) 










result 


ae 


llUg^, sjnjx; 


m 




le 


^^■fej'- result; 


ti- 
BS 


B^^OTp = sin X * cos.h + cos_x'-^^^g 


m 


cosj = cos.x * cos.h - soji.x * "^^m 


!t- 








3) 











you think? We haven't really changed 
the nature of the series or even saved 
any storage or arithmetic requirements, 
but the series certainly looks simpler, 
and the malfa is easier to woric widi, 
assuming your programming language 
supports complex arithmetic. C does 
not, but most C++ vendors are supply- 
ing a complex number class as part of 
their libraries, and ANSI is working to 
establish a standard class. See P.J. 
Plauger's "Complex Made Simple" for 
the latest information.^ 

At flist glance, it may seem that 
we've required even more storage 



because the index of Cn must nm over 
all mtegers, where that of a^ and b„ ran 
only over the noimegative ones. 
However, the c^ are not all indepen- 
deet TimvsJiaes for positive and nega- 
tive n aie ewples ccmji^ates of each 
other: 

e_, « e,' (15) 

So ^ used te evaluate and store coef- 
ficients only for nonnegative values. 

Let's derive expressions for the 
coefBcia^ Cg. We can get these 
expressions directly from the defini- 
tion of the series, as we originally 
fisand aad ^ tot yen might be more 
convinced we're on the right track if 
you see them derived from Equations 
11. For Co, we have: 



Co = ^ / fix)dx 



EquatitM 5 gives tie 

a^. Substituting the ii 
Equation 8 gives: 



-*-jr 

a, = ^ / fix){e^ + e-^)dx 
-It 

Similarly, for hj,: 

+11 



for 
from 



Now, coc 
we get: 



i ttese » Iptdon 11, 



Cn = 



J_ 
2.T 



-n 

ii-i) I f{x)(e^ - e-^)dx 



(Watch out for all those -i terms. Their 
product is minus one.) Combining 
terms gives: 
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e^dx 



L 



J 2f{x)e-^ 



dx 



(17) 

Similarly^ the case for n < gives: 

which becomes the same as Equation 
1 8 if we replace -n by n. Since this for- 
mula also matches the special case n = 
0, we end iq> with cms expxes^tm for 
all n. 

We now have a rigorous definition 
of the Fourier series in its conq)lex 
form. It is: 



fix) = Y^c„e" 



(18) 



(19) 



Armed with suitable classes for 
complex numbers and their vectors, we 
can rewrite the code of Listings 1 and 
2 to use these equations. The results 
are shown in Listuig 3. Small, aren't 
they? Fvmction fourier now contains 
eight executable lines of code. That's a 
lot of processing in a small package. 

IS IT DFT YET? 

The listings deserve a bit of 
explanation. First, I've clearly 
defined a class C.Vector that is 
virtually identical to class Vector, 
except its elements are conq)lex. I'm 
using the Borland class for complex 



USTING3 

HFFT. complex version. 

fci// Evaluate Fouriei^^ifis^ 
| doubiLe<fourier(doufi]je''3(,%;.yector k c){ 
I?; intm = c.rankO; 
f conplM-JW '= «xp(con|iLKt(0, x)) ; 
^ coflifilex wi = conpljex(l;iO}j^. , . 
J double result =|reail(c[(fl)i"/';^^;''''' 

-forCjjtt i=l; a^ 
I* "wi *= «; 
f result += 2.0 * 

L return result; 

\~ 

h/l CMipute Fourier coeffocients? 
^-void fourfit(Vector ftf , C.Vector kz){ 
If.rankO; 
.rankOj.;,,:^ 

TiV = exp(coB^x(0,\2.C 

yi = confiLextt^l 

!X «1| 

r(uit.i=0; 1 < n; 

pj =:COfflalsx{l; 0);' 
|or{int j=0; i < m; j^-^){ 
|c[j]^=f[.l > 




// Regenerate Fourier function 
void ajiv_four(CJector k, Vector tf){ 
int n = f.rankO;- 

int n = c.rankO; ■ 



.2.0 * 





complex la = 

cofnplex yij; 
f = 0; 

fordnt 1=0; i<n; i-h-){ - 
wij - coinplex{l,.;0)j:V 
for(int j=0; 'J<l^;^g^-^){ 

f[i] += 2.0 * real(c[j] * idj); 

yij *= ui; 

Hi *= 



} 



mm 



numbers. In fourier, I generate the 
fimction value on the assimiption that 
it's real, which was the same assump- 
tion Fourier made. Because of the sym- 
metry in Cn for this case, we don't need 
to store the elements for negative n, 
which saves us half of the storage. 
Finally, take a look at the listings for 
fourfit and inv.four side by side. 
You'll see that they're virtually identi- 
cal. So, do we finally have a discrete 
Fourier transform? Yes and no. In one 
sense, we have more than the DFT 
since these functirais are more flexible. 
Remember, they were originally writ- 
ten to support the concepts of the 
Fourier series as Fourier himself 
defined it. That is, we express a peri- 
odic fimction f(x) in terms of its sinu- 
soidal components, with the purpose of 
approximating the original fimction or 
otherwise doing tricks with it. 

Remember, fimction fourfit was 
designed to be a many-to-few transfor- 
mation; I expected to generate a set of 
coefficients much smaller than the 
nimiber of points in the original fimc- 
tion. Function fourier was a few-to- 
one conversion, but there's nothing 
except Shaimon's limit to keep us fi-om 
using the routines with larger vectors. 
That limit requires that the order of c 
be no larger than half the nimiber of 
samples. Since c is complex, the nimi- 
ber of elemmts in c is precisely the 
sams as those in the original fimction 
vector. A typical DFT algorithm is 
written based on the assumption that 
this will always be the intended use. If 
this is what you want, and you're only 
interested in dealing with real func- 
tions, you can simply declare vectors f 
and c with dimensions n and o^, 
respectively, and you've got your 
DFT. 

Real fimctions, did you say? Are 
there any other kmd? Of course there 
are. Although we've been working 
with f(x) as though it were always real, 
there's really nothing in the math that 
depends on that assumption. Equations 
18 and 19 work just as well if f(x) is 
complex. The only problem is that the 
elements of c for positive and negative 
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If we can't yet call 
the funcfiens I've 
shown DFTs, It's 
because we 
haven't yet 
defined the Fourier 
tM^orai Hsdf. 



indices will no longer be complex con- 
jugates, so we must store the values for 
the full range. This, in turn, gets us into 
the issue of how to organize the vector; 
C++ doesn't allow negative indices. 
Programming issues aside, but ttte 
approach is practical. Although we'll 
now have twice as many elements to 
store, we'll still have an equal number 
in both vectors. Many DFTs are writ- 
ten to work with complex functions. 
The only problem witii that approach is 
that 99% of the time, you'll be working 
with functions that are real, so you'll 
be wasting half the data storage (and 
using four times the computing time) 
for no good reason. This causes a cer- 
tain bother. 

I've seen programs that solve the 
dilemma by transforming two real 
functions together, putting them in the 
real and imaginary parts of the com- 
plex vector f. But that approach strikes 
me as something of a kludge — a way to 
do something with the storage that's 
sitting there, allocated. How often can 
I expect to be dealing with two totally 
different functions, but two I happen to 
have tabulated with the same number 
of data points? Not very often, I'll 
wager. For most cases, I prefer the rou- 
tines of Listing 3, which assume that 
f(x) is real. However, the symmetry 
between fourfit and inv.four is more 
nearly perfect if we use complex val- 
ues for both vectors. 



HOME ON THE RANGE 

If there's a reason we can't yet call 
the functtons I've shown DFTs, 
it's because we haven't yet defined 
the Fourier transform itself. It's about 
time we did. 

To do so, we must first talk about 
the range over which the fimction f(x) 
is defined. As we've been using the 
Fourier series, we've taken that range 
tobe-n .. +71. Such a range is reason- 
able yabm we're talking about an angle 
that's naturally repetitive over that 
interval. Even then, we might prefa: a 
range beginning at rather than centered 
on zero. Again, in the case of an angle, 
it really makes Uttle difference which 
range we ehoose, since the function 
repeats itself on an interval of 27r in 
any case. As a start, suppose that the 
fimction given to us is not centered at 
the origin, but begins there. We can 
easily translate the range by a substitu- 
tion of variables. Let: 

y = x + ;r (20) 

so, equivalently: 

x = y-;r (21) 

Let's define a new function, g(y), such 
that: 

g(y) = f(y-jr) (22) 
Our Fourier series b%tMU»e: 

+00 

giy) = f{y-n) = Y^c„e''<y-''y 

The last exponential is independent of 
y, so we can tliink of it as merely mul- 
tiplying c„. Define a new coefficient: 




c„ = c„e 



(23) 
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Then our series looks exactly like il did 

before, except for a prime over the 
coefficient, as follows: 



BOSTON 

SYSTEMS 

OFFICE 

— m 1 

TASKING 

Fy6upfreesdein<Ktdislt ard^ 
our competitive B'erfcBmarlc^l 
1-800-458-8276 
Fax 617-320-9212 ' 



TASKING Europe -1-31 33 558584 
TASKINCJapan +81 334 050511 
M tndennkBd naine are the propeity of the respective ownen 



CmCUE #24 ON REABra SERVICE CARD 



+00 



giy) = 2^c'„e«y 



(24) 



That prime is a hint that the coeffi- 
cieat anay will be different. The dif- 
fermce is embodied in Equation 23. 
From the definition of c-, we can write: 



Cn = 



or. 



'dx e 



Maicing the same substitution for y and 
W^ig, toe values of ywbim*. Is at its 
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(25) 



This equation also has precisely the 
same form as Equation 19, except for 
tbe range of integration. The implica> 
tioaa isthis: The algorithm doesn't real- 
ly care what you define the range to be 
(even though a range is implicit « lite 
formulation of the computer pro- 
grams). After all, a periodic function is 
a poiodic function, meanii^ that it 
will go on forever in both directions, 
repeating at intervals of 27r. The point 
you pick to be the starting point is 
rather arbitrary and depends more on 
your convenience than any deep math- 
ematical requirement. Just be consis- 
tent, and everything will be fine. 

On the off chance that you do need 
to change the coefficients after they 
were generated. Equation 23 gives 
what you need. In general, an offset 
throu^ any an^e is equivalent to a 
rotation in the complex plane. But a 
rotation through an angle like n is very 
special; it really isn't so much a rota- 
tion as a simple sign change. This is 
most obvious if we look at the real 
rather than complex series. We know 



sin<x -i- ym'ya^ite JEM rm-Hm xsi&m 



Similariy: 

eos(jc + tm) = cosxcosiiff 

Finally, cos wt obeys the simple rules: 



cos mr = 1 
cos mr = -1 



(n even) 
(n odd) 



To get one set of coefficients from 
the other, change the sign of every odd 
^Msffident We em mt iiat tiiis 
scheme works for the special case of a 
square wave. If we define the function 
as odd (by p^tt^ ibe tayi^ttet at the 
origin), we get no cosine terms, and the 
sine terms are zero except for odd 
orders. Thus, every ncmz^ coefficient 
of the square wave gets negated when 
we shift the origin by tt. And since 
evey term is n^ated^ tta@ resulting 
value of the function as evaluated 
using the series is negated as well. 
Hiat's precisely tiie result we ^ould 
get, since shifting a square wave by t 
is indistinguishable from simply negat- 
ing it. 

UW OF AVERAGES 

Thapot^^^ lUs di^ussion, one 
point was probably painfully 
obvious to the mathematicians 
ttaoag m, \k& p«faipl mm so for the 
rest. Going back t© (fee Heal series once 
more, consider the d^oition of ao, or 
rather, <^2- ^^^ota Equatim S: 



ao=iffix)dx 



so: 



(26) 



Here we're mtegrating over a range of 
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In, then dividing by that same factor. 
Mathematicians will recognize imme- 
dii^ly that diis integral is the simple 

average of f(x). So the constant term in 
the Fourier series is simply the average 
of the function. Similarly, every other 
coefficient is also an average. The for- 
mula for the complex coefficient: 



also represents an average, this time of 
fteimction: 

This inno- nature of the integrals 
will be important when we deal with 
intCTvals other than In, which we'll do 



FUNCTIONS OF TIME 

The range -n..+it may be all 
right for angles that naturally 
repeat on that interval. But a 
goodly portion of our work ends up 

dealing with functions of time, and few 
are going to be so accommodating as to 
repeat every 6.28 seconds. To deal 
with functions of time (or any other 
measurements than angles), we need a 
much more flexible definition of Ae 
Fourier series and coefficients. 

Fortunately, we can easily accom- 
modate the changes based on what 
we've already learned. Consider a 
function of time with period T. Let: 



x = (o t 



(27) 



and we'll find out what value w must 
have. For things to work properly, we 
require that when t = ^/2, x will be n 
(and when t = T, x will be 2jr). Solving 
for CD gives us: 



(28) 



The ratio Vj is called the frequency 
f, and many formulations of the 
Fourier transform use f instead of o). 
I'll use f sparingly here, partially to 
avoid concision between tbe fi«quency 



f and the function f(x). Substituting our 
new definition for x into Equations 18 
and 19 gives: 



(29) 



^^lece: 



m 



Looking at Equation 29, we see that 

each term is a sinusoid with a frequen- 
cy of n times the fundamental frequen- 
cy. Each coefficient c„ defines the 
amplitude and phase of the harmonic 
associated with its corresponding fre- 
qntm^ nf. In other words, gives us 
the spectrum of f(t). 

With that thought, we begin to see 
the concept and value of Fourier analy- 
sis. This is no longer about approxi- 
mating a fimction. Founer analysis is 
about analyzing a sipial to ^imnine 
its frequency content. That can be 
important in a wide variety of fields, 
from communications to stereo to 
image processing to data compression. 

For these kinds of applications, the 
set of coefficients c„ become fer more 
than just the values we use to regener- 
ate the signal. In modem analjrsis, we 
often don't bother to regenerate it since 
we already have the original. For spec- 
tral analysis, the coefficients are the 
message, tat thing we'le talasiled in 
seeing. 

The change from using the Fourier 
series to doing Fourier analysis is both 
trivial and profound. It involves a shift, 
not so much in how we do things, but 
in how we think of the results. When 
we're using the equations to approxi- 
mate a function using the Fourier 
series, the vector c is siiiq)ly a set of 
coefficients whose values we need to 
get the right fit. When we're doing 
Fourier analysis, c is the desired out- 
put or spectrum that describes the 
function. 

Look at it Ms way — we began with 



a continuous fianction of time, f(t). We 
expressed that function as a table of 
values or, more accurately, a vector. 
Strictly speaking, once we've convert- 
ed the fimction to tabular form, we can 
no longer think of it as continuous, but 
we do. On the other hand, our algo- 
rithm produces the output vector c, 
which our brains still want to think of 
as being a discrete set of coefficients, 
mere byproducts of the algorithm. But 
just as f(t) can be a continuous func- 
tion, a table of values, or a vector, we 
can think of c as being a fimction in its 
own right — & fimction of <w rather than 
time. Just as we've seen the Laplace 
transform convert a fimction from the 
time domain to the s-domain or the z- 
transform to the z-domain, we now 
have another domain to play with: the 
frequency domain of Fourier analjrsis. 

Over the years, people have come to 
adopt the imiversal symbol h (or H) to 
represent transforms. We can think of 
the transformed fiinctions not so much 
as new functions, but as different rep- 
resentations of the same function. 
CiHic^&ially: 

h(t) <*> H(w) H(s) <=> H(z) 

Which representation we get depends 
on which transfonn we applied to get 
it. If nothing else, this notation, in 
which the nature of the function is 
implied by its argumrait, saves time 
trying to think up new symbols for 
each representation. 

OK, so we can think of c as beii^ a 
function of w. Does that mean we can 
also think of it as being continuous, 
like f(t)? Well, no, we can't, not when 
it's derived from a periodic function. 
We have no need to do so, since the 
only frequencies that appear in the 
function, appear at intervals of the fiin- 
damental, which we might call a>o. In 
the jargon of spectral analysis, there's 
no energy at other frequencies, so the 
function representing the spectrum 
really is discontinuous; energy ^pe^ 
as infinitesimally wide spikes at the 
harmonic frequencies, and the wapW- 
tude is zero, otitasrwise. I%ure 1 shows 
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a plot of the amplitude vs. frequency 
for the square wave we've used as a 
test case. In Figure 1, I've normalized 
the frequency to the fundamental a>o. 

Still, Hus whole new way of looking 
at things certainly raises our curiosity: 
Do ftmctions exist for which H(a)) is 
continuous? The aiiswer is, yes indeed, 
and that brings us to the notion of the 
Fourier integral. Not surprisingly, if 
periodic functions lead only to func- 
tions in the frequency domain that are 
discontinuous, we must look to the 
realm of nonperiodic fimctions for 
mything else. It's interesting, though, 
that functions exist that have only dis- 
crete frequencies, but are still i^ioiod- 
ic. Hamming gives a very m^gki 
example:^ 



/(/) = sin<a t + sia-^joTt 



(31) 



As you can see, iliis fiuictioD em- 

tains only two pure sine waves. In a 
plot of the frequency spectrum, it 
would have only two lines where the 
values were anything but zero. But the 
two frequencies are not commensurate. 
Therefore, once started up with a given 
set of initial conditions, the ftmction 
will never, ever, repeat itself Each 
time one of the sine waves passes 
through a certain phase angle, say zero, 

FIGURE 1 

Spectrum for a square wave. 



the other wave will have a different 
phase than it has ever had before. This 
ftmction is not periodic. 

We can never represent such a fimc- 
tioQ as a Fourier series because no 
combination of harmonics can ever 
represent both waves. If we try, we'll 
end up with a discontinuity at the edges 
of any interval we try to define. And as 
you know, discontinuities lead to ring- 
ing, Gibb's phenomenon, and so on. A 
Fourier analysis, based on any finite 
interval, will yield frequency compo- 
nents that doB't im&f exist in tlie orig- 
inal ftmction. 

If we're to have any hope of repre- 
sentii^ a fu^^iitt Mlas tbis one, we 
must turn to a representation that 
includes all frequencies, not just those 
that are haimonics of a fundamental. 
And that can only be done with an mte- 
gral, not a summation. 

THE FOURIER INTEGRAL 

The question is, how do we get 
from a series to an integral rep- 
resentation? Does such a thing 
exist, or are we just dreaming? I can 
give you four approaches. Each of 
them has some merit, but none of them 
completely satisfies. Together, they 
make a powerftil argument. 
The first argument recognizes diat 
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sine waves represent the solution to the 
differential equations describing a sys- 
tem. The principle of superposition 
that holds for any linear system says 
we can mix the solutions to get a gen- 
end solution. Take a vibrating string, 
for example. Typically, the general 
solution of that problem results in a 
variety of sine waves, all rippling back 
and forth together without interacting. 
The system behaves as though each 
sine wave were doing its tiling inde- 
pendently of the others. Now, it might 
seem that a vibrating string is a poor 
racample since we all know that tiie 
solutions for such a string lead to only 
certain frequencies, a fimdamental and 
its haimonics. But Ibat's only true fat 
the unforced motion. We could hook 
up a piezoelectric transducer, for 
example, and drive the string at mxy 
frequency we can dial in on an oscilla- 
tor. Indeed, the equations of motion of 
the string don't preclude sine waves of 
any frequency — only the boundary 
conditions force the harmonic solu- 
ticms, and that's typical of a dynamic 
system. 

The second argiunent is similar. 
Suppose we have some medium, iay a 

piece of Styrofoam, and we want to 
study its sound transmission qualities. 
Better yet, make it a telephone line. 
Physically, we know tiiat we can hook 
up a transducer at one ead and a micro- 
phone at die aOier and drive tiie trans- 
ducer with any frequency we dial in. 
We can plot fbs frequency response of 
Uiis medium by sampling its re^ionse 
at different frequencies. Unless it's a 
very special medium indeed, we know 
deep in our l^urts that the frequency 
response curve is going to be continu- 
ous, not a series of lines like we see in 
Figure 1. If tiie physical system can 
transmit frequencies of any value, 
surely the mathematics must support 
them as well. 

The neatest approach for justifying 
&e integral comes from Hamming. 
The general idea is to ^proximate an 
aperiodic fiinction that we don't yet 
know how to analyze by using a peri- 
odic one diat we do. We'll see what 



happens as the approximation gets bet- 
ter. Consider, for example, the unit 
pidiB ■^bmm. it Figure 2a that has a 
poise milk et t. ify&maticany: 

f(t) = l .«2<t<V2 

f(t) = othermse (32) 

As a flnt mi very rou^ approxima- 
tion to our pulse, we'll use our old 
friend the square wave, appropriately 
o£^ as m MgWe 2b. This Amotion 
matches f(t) perfectly, for that one lone 
pulse, but of course it repeats else- 
whore cm tiie tee 1^ a period T 
equal to 2r. We'M ^ this fimction 
f2(t). All in all, it's a psily bad fit to 
flit), but it's a start. 

As the next step, let's make the 
function longer by sliding out the adja- 
cent peaks on the wave, as in Figure 
2c. This is not equivalent to just scal- 
ing the wave. We're not changing the 
pulse width, only the time betwerai 
pulses. The period increases to, say, 4t. 
We'll call this function f4(t). 
Eventually, we'll end tap wift a whole 
set of such functions, each representing 
better and better approximations to f(t). 
The phantom pulses demanded by our 
condition of periodicity go sliding off 
the page in opposite directions, head- 
ing for infinity, lotving oiu' original 
pulse standing alone in all its pristine 
glory. Somewhere out there, those 
pulses still extst-^an Mnite number 
of them, in fact. But for practical pur- 
poses, they end up so far out they 
might as well be gone. I'm going to 
assert that in the limit they are gone, 
and that limit occurs when the period 
f^i isat approxnnating fym^m tesshes 
infinity. 

The fourth and final justification is a 
mathematical ^e and should be 
enough to convince even the most 
hardened skeptic. We'll begin by sub- 
[ 30 into Equation 29: 
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FIGURE 2 

Approximations to unit pulse. 
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Next, let's relabel o) to (Uq, and recall 
itidefimtion: 

Also, we'U define eo to be: 



(34) 



can think of &e c(»»liaiit oio as lidng a 
Ao), the step between summation 
points. This is a crucial step because it 
ghnsB us stnned^^ a ^p t ea dfeii^ an 
infinitesimal. As T grows larger and 
laigCT, (Oq = Aeo must necessarily grow 
smaller. Moving things anmA & bit 
gives us: 



The reason for this last step is that 
we want a variable other than n that 
changes for the different terms of the 
summation, coq is a constant (for a 
given T), but o) is clearly not. 

With these definitions, Equation 33 
becomes: 



f(t)=j2 



-5^ 



(36) 



J 'dt 



(35) 



Next, take a look at the last expo- 
nential. This is changing with each 
increasing value of n because of the 
increase in co. And how much does <u 
increase between successive sums? 
Why, it's by a>o (ree Equation 34). We 



Now we have something that's truly 
beginning to look like an integral. 
We've managed to eliminate the inte- 
ger n everywhere but in the summation 
itself Our multiplier looks an awful lot 
like an infinitesimal at the end. The 
next step is to consider what happens 
to the formula as T increases. A<w must 
go to zero since it's proportional to V^-, 
giving us two ingredients we need to 
proceed to the integral: A decreasing 
infinitesimal, and points that are 
becoming rdate closely packed. Hie 



only thing that can stop us now is if the 
range of n somehow gets limited. But 
that's far fixan the case because not 
only is n proportional to w (see 
Equation 34), but the proportionality 
constant is T: 



CO _Tw 



(37) 



While w passes through its range, n 
goes faster and faster to its infinities as 
T increases. Not only are we still going 
to get a sum over all values, we now 
know that the new viffiable of integra- 
tion is CO. 

In the limit as T approaches infmity, 
our sum becon^s an integral: 



(38) 



Now we can recognize the two trans- 
forms, so we can pull them apart again. 
We end vap with a formal definition far 
the Fouria trai^oim ^d its invra^: 



hit) = j H{(o)e^'dw 

eo 



(39) 



(40) 



You might have noticed that the place 
^ere we put the factor of 2n was pret- 
ty aibitiary. We could just have easily 
have put it in the first integral, and 
some people do. We could even put it 
both places (using its square root) to. 
get the ultimate in symmetry. But if we 
want to talk about fi-equencies and 
spectra, it must go where I've shown it. 

SERIES OR TRANSFORM? 

Like one of those pictures that 
look like different things 
depending on how you look at 
it, we can think of the two integrals in 
Equations 39 and 40 in two profoundly 
different ways. On one hand, Equation 
39 is merely an extension of tiie 
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Fourier series and lets us approximate 
a function (not necessarily periodic). In 
lluit worldview, tibe companicm equa- 
tion, Equation 40, is just the way we 
iind the coef^cient function. 

On the oflier hand, we can think of 
the equations as defining a transform 
pair just like the Laplace transform. (In 
fact, the Fourier transform came first.) 
Just like the Laplace transform, we can 
use the Fourier transform in all kinds 
of interesting ways. In this wtwldview, 
the transforms give us the mechanism 
for moving from the time to the fre- 
quency domains. Operating on the 
assumption that two domains are better 
than one, and three or four are better 
than two, tfiis certainly seems a plus. 
What might be difficult in one domain 
may turn out to be trivial in another. 
The more we have access to, the es^tse 
life is likely to be. 

Ironically, now that we've gone to 
aH &at trouble to justify the defmition 
of the Fourier integral and Fourier 
transform using continuous fimctions, 
we ffiust now go back the other way 
because we still can't represent the 
fimctions in a computer except as 
tables of values. In oflier words, we're 
going to end up right back where we 
started, with the routines of Listing 3. 
But we now know some important 
things that we didn't know before: The 
procediures can be used for nonperiod- 
ic functions, and the frequency spec- 
trum is really a continuous function, of 
which we're seeing only a part. We'll 
continue in that vein next month. 
Before we leave the world of continu- 
ous fimctions completely, you may be 
wondering whatevo' Yisppaaed to our 
imit pulse. We now have the tools to 
compute its transform, so I'll leave that 
wift you as a parting gift 

Recall, Hie fimctioQ was defined to 
be: 

f(t)=l -V2<t<V2 
f(t) = otherwise 

Inserting this into our new defmition of 
the transform reduces the range of the 
integral to be: 



HGURE 3 

Transform of a unit pulse. 
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-'A 

The ta^pal is a fundamatt^ ome: 

'A 



(41) 



■ f iwz UoT \ 



But this is ahnost i 
nition of the sine fm/^^BM 

sin(x) = ^ 

So our result is: 



je-^'dt = lJ^^ (43) 



-y2 

Our transform, then, is 
F(fl,) = ^smfc 



(44) 



Because of the presence of (v in the 
denominator, this function has the 
form of the familiar sine integral: 



smjc 

X 



(45) 



A ffSB^ of tile transform is plotted in 
Figure 3. 

That's all we have room for this 
month. We'll wr^ up this series next 
month with a look at tibe fast Fourio' 
tnmsfonn. 

Jack W. Crenshaw is a sU^ scientist at 

Invivo Research in Orlando, FL. He 
did much early work in the space pro- 
gram and has developed numerous 
analysis and real-time programs and 
holds a PhD in physics from Auburn 
University. Crenshaw enjoys contact 
and can be reached via e-mail at 
72325.132 7@compuserve. com. 
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AU About 
Fourier Analysis, 
Part III 

The series continues with the theory 
and software foundation you'll need 
to understand and implement Fourier 
transforms. 



If you've been following this 
series, you have noticed 
that we have seemed to oscil- 
late between theory and coding. 
We began with a very mathe- 
matical and theoretical treatment of the 
Fourier series and then swung the 
opposite way to siis^le and practical 
code implementations. 

Last month, we continued in the 
software mode and discussed ways to 
generalize the previous implementa- 
tions before returning to a very mathe- 
matical treattti«»t, including the moti- 
vation for and derivation of the Fourier 
transform. Here, we'll cover software 
implementations ifee Fourier trans- 
form. But first, we'll stick with the 
pure mathematics and play with the 
notions of the Womm tmosfonn in its 
theoretical guise in order to get a good 
feel for how it works and what it's 
good ibr. Next month, we'll wrap the 
series up with a look at the all-impor- 
tant fast Fourier transform (FFT) and 
its implementations for both real and 
complex functions. Hang in there — 
we're almost done. If everything goes 
well, you will come out of this series 
with a set of software implementa- 
tions in your pocket and. more impor- 
tantly, the knowledge of how and 
when to use ttem. ( 



At the end of last month's article, 
had just introduced the formal defini 
tion of the Fourier trm^^flm §m: 

hit) = / H{a>)^do} (1 

The symmetry between these two in 

grals should be obvious. Except for 
sign and the factor of 2;r, the forms 
identical. In fact, most compu 
implementations use the same code 
both transforms, with only some p 
and post-adjustment for the inve 
transform. Once again, this transfo 
represents the logical extension of 
Fourier series, in complex form: 

fif) = 2^c,e^ 

-Vi 

Remember, the Fourier series 
viMlcs for imctions that are peri 



r 



with a finite period T. Last month, we 
saw how the Fourier transform arises 
as we take the limit as this period goes 
to infinity, thus allowii^ us to handle 
nonperiodic functions. 

To get a feel for the way the trans- 
form behaves, let's look at a few prac- 
tical examples. The most simple func- 
tion of time that comes to mind is a 
craistatit: 



h{t) = 1 



(5) 



You probably know what this function 
would look like represented as a 
Fourier series (not transform). All 
coefficients vanish except for the one 
constant term cq, so the tx^foim of 



(6) 



But what does it look like when we 
apply the rules for the integral trans- 
form? Substituting the constant for h(t) 
into Equation 2 gives us: 



(7) 



I see two problems with this inte- 
gral. First of all, for any nonzero value 
of CO, the integrand is simply a sine 
wave at that frequency. Now, we all 
know that a sine wave oscillates about 
zero, and its average value averaged 
over a fidl cycle is zero. And an aver- 
age is exactly what we have for the 
finite series coefficients defined by 
Equation 4. But in this case, we're not 
exactly averaging, since we aren't 
dividing by the range. VIore impor- 
tantly, our range of integration is far 
from a period; it's infinite. I could 
wave my arms and assert that the for- 
mula is supposed to mean an infinite 
number of full periods, but I suspect 
this argument will leave you as unsat- 
isfied as it does me. .A.fter all, that 
range is different for each frequency. 
It's extremely difficult to argue that 
we should use different "values" of 
infinity for each frequency. 



You will come out 
of msMiiSwIih 
a set of software 
implementations 
to iour pocii^ ami 
the knowledge of 
km and when to 
use them. 



The ^ond problem is vfbait hap- 
pens for <B = 0. At that fi%quency, we 
get: 



(8) 



It would appear that our o^sform 
of a simple constant function doesn't 
yield a nice constant after all. but 
rather a value of infinity at the origin 
and zero elsewhere. And even that zero 
seems a little suspect because of the 
need to assume that the integral of all 
sine waves is zero. So our simple little 
first example is not so simple after all. 
When we understand what's going on 
for this case, we'll have a much better 
feel for wiat the Founer transform is 

irSALLMTIIEOIirrS 

We can get a first hint as to 
the meaning of our result by 
looking at the units of each 
fiinction. Of course, the function h(t) 
can be used to represent any variable 
that's a function of time, including 
voltage, current, distance, velocity, or 
anything. In the case of the Fourier 
series, the sine-wave multipliers of 
each coefficient arc Jimensionless. so 
whatever units /in has, all the coeffi- 
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dents c„ have the same units. 
But look at the integral in Equation 

2. The factor e"»" is dimensionless, but 
our function h{t) is multiplied by a time 
variable dt. Thus, wbatever units h(t) 
has, the units of H((o) camiot possibly 
be the same. If h(t) is measured in 
volts, H((o) must have units of volt- 
seconds. Or, since the units of w are 
time-i, we can think of H{w) as repre- 
senting volts per unit frequency. (Here, 
I'm using the term "frequency" a little 
loosely. Strictly speaking, co is not fre- 
quency but 2jt tim^ die frequency. M 
casual conversation, we tend to use the 
two terms interchangeably.) 

B(o)) represents a density func- 
tion — the amount of energy contained 
in the signal having a frequency 
between co and ort-efco. Faf menit sig- 
nals, certainly those that are truly non- 
periodic, this amount of energy will be 
proportional to da>. But tiiat assump- 
tion breaks down when we have a truly 
periodic signal. In that case, all the 
«iergy occurs precisely at spedfic fr^ 
quencies. The more we try to narrow 
down the precise frequency by reduc- 
ing dco, the greater the density s^peais 
to be. In the limit, as dco goes to zwo, 
the density is indeed infinite. 

We get a similar result if we com- 
pare the graphs of the transform for a 
periodic and nonperiodic fimction. 
Last month, I showed the graph of the 
Fourier coefficients for a square wave, 
reproduced here as Figure 1. The 
square wave is periodic, so we must 
use the Fourier series. We get a collec- 
tion of spites at specific frequencies, 
as you can see. (In reality, die eoe^ 
cients of even and odd frequencies 
alternate in sign but we usually plot 
only the absolute value of the coe£5- 
cients since, in gmeral, Ihey may be 
complex.) 

Last month, we also generated the 
Fourier transform (using the integral) 
for the nonperiodic function represent- 
ing a single unit pulse. We found this 
transform to go like Ac sine integral: 

sin.r 



the absolute value of which is repro- 
duced Figure 2. Comparing die two 

graphs, you can see that they have a lot 
in common. In both cases, the ampli- 
tude changes as i/o). But there is also a 
profoimd difference. In Figure 1, we 
have clean spikes of finite height but 
zero width. M Figure 2, the spikes are 
smeared out into a spectrum. Each lobe 
of Figure 2 has not only a finite height 



but also a finite integral, as represented 
by Ae area imder die curve. But the 
spikes of Figure 1 have no "area under 
the curve" at all — their integrals would 
be zero since they have no width. This 
just serves to emphasize that the 
graphs, though they might look the 
same, aie really plots of totally differ- 
ent things. Figure 1 is a graph of ampli- 
mdes, while Figure 2 is a graph of fre- 



F1GURE 1 

Spectrum for square wave. 
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Ihmsform of unit pulse. 
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quency 4ssslty. Since die ^ikes of 

Figure 1 are infinitely dense, we 
should not be surprised to see them go 
off ^ seale when plotted the same 
way as in Figwe 2. 

iHE BRiiG DGLm Fumnioii 

All of this will become clearer 
with the introduction of the 
Mme (teita function, always 
written i5(x). This peculiar function is 
zero when x is nonzero and infinite 
v^m, X m mm. Bat it is ixuBmte in a 
speckd wxf, so &at its a^BB§0M* 



1 



(9) 



You might have noticed Itat I #ia't 
^ow any oa the w6m§sA, 'Qist's 
because it doesn't matter what ffei^ pi 
as long as the range includes ttomte 
of zero, tie iol^pil dBi mubees 
is no different than, say, the integral 
from -0.001 to +0.001. 

Becsise tMs starage fbt^tiba is 
nonzero only for one specific value of 
X, it serves as a filter or selector for 
anything else in im^piL HffiK 



fsix)f{x)dx = m 



(10) 



In controls jargon, the delta fimction 
k a wait knpulse. It occurs over zero 
time, so it has infinite amplimde but a 
finite total effect because its integral is 
finite, in fact unity. Such impulses 
can't really exist in nature, but we can 
come very close. For example, when 
two billiard balls collide, the force 
imparted on them is very large for a 
very short length of time. To compute 
^ motion of axe billiard balls, we 
iM*t to know the details of how 
tMs fcxax changes over time. The time 
the balls are in contact is so short that 
they don't really have time to move 
much. For all practical piuposes, the 
(iaage in velocity (and spin, if you're 
a pool shark) of each ball occurs 
instantaneously. The unit impulse is a 
mythical but useful idealization of 
real-world events. 

Mathematically, the delta fimction 



can be thought of as die derr^ve of 

the unit step shown in Figure 2. In the 
unit step, the input (force, voltage, or 
whatever) is zero M affl Wm up t» tte 
point t=0. At that time, it changes by 
one unit, whatever the vmits of the 
input ate. Clearly, ^ d»i«^ve of this 
input, which is the slope of the curve, 
is infinite at / = 0. And yet, equally 
cleffitff , iti mt^ral h finite «ice %y 
definition the integral of a derivative 
yields the original fimction. Thus, the 
kipolse is precisely die <krivative 
of the imit step. 

Let's see what happens when our 
IB^ teetiEHi is a (^tta tacticm (nr unit 
iaq>ul8e. From Equati<»a 2, we have: 



H(.eo) 



1 

2;r 
1 

" 2jr 
or simply: 



/ S it)e- 



"dt 



(11) 



in odier words, when our input is die 

unit impulse, the fi-equency spectrum is 
a constant. This is another way of say- 
ing that the unit impulse contains all 
frequencies from D.C. to gamma rays 
with equal amplimde. 

Wie get an even more interesting 
result when we perform the inverse 
transformation. If our fimction in the 
fiiequency domain is a delta fimction: 



H(a>) = & (co) 



(12) 



(implying, of course, a single fireqiten- 
cy of zero), we get: 



^<f)=l 



(11) 



which is precisely the cra^ant value 

we had trouble with going the other 
way. While a constant value may give 
us preplans in die mtegration over 
infinity, the inverse transform is trivial. 
Fortunately, we have the formula for 
both transforms, so we can go in 

whichever direction is the easiest. 
What about a single frequency? That 
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one's easy too, because we can trans- 
late the delta function to any point 
other than the origin by a simple 
(tenge of variables, lliere's only om 
small problem: In the frequency 
domain, if we have a spike at some 
<9alue, say Q, we must have a ccme- 
sponding one at the negative frequency 
as well. So our fimction must look like 



MO = cos£2f 



(15) 



H{q)) = ^[S ia)-a) + S{co + Q.)] 

(14) 

Taking the inverse transform gives us 

— oo 




(You might be wondering why we got 
a {KDce eeskne wave with no sine wave 
or phase shift. That's because of our 
choice of multiplying constant, unity. 
M mgi ap^ki^ims, we'd multiply die 
amplitude by a constant which might 
be complex, leading to an arbitrary 
^ase im^e.) 

WHY TRANSFORM AT AU? 

At penM, stay feel that 
you've learned more than you 
care to know about Fourier 
integrals. But one Mm wt sM havea't 
discussed to any degree is a very fun- 
damental question, what the heck are 
they good for? How do diey help me 
get my problem solved? 

A hint of the answer comes irom 
considering how other tnmsforms are 
used. We discussed &e use of L^lace 



transforms in my article. The A^^ft- 
bet fi-om S to Z." (Aug. 1994, pp. 60- 
73.) The idea goes like diis: Suppose 
we have a dynamic S3rstm dssoaM 
by a linear differential equ^oa: 



„ d''x d''-^x 

OoX = 11(f) 



(16) 



where u{t) is a known but arbitrary 
forcing fimction. We'd like I© knew 
the behavior of x as a fiinction of time.. 
The key word in the previous sentenee- 
is linear, mi that point bears empha^ 
Everything we've done here, including 
the Laplace and z-transforms, is based 
on the assumption that the different' 
equation is linear, meaning that x and 
its derivatives enter only as additive 
trams. For linear vepgekmOk ^ |nin^>j 
pie of superposition ^li<a^ ^ a fiinc«| 
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tion jc(0 is a solution to the homoge- 
neous equation (where u(t) is set to 
zero), then 2x(t), 3x(t), or any constant 
Utem x(f) is also a solutiaa. If any two 
functions are solutions, so is iiie siHn of 
those functions. 

Transfonns like the Fourier trans- 
form and Laplace transform are useful 
to us only when the differential equa- 
Hon is lin^. In fact, for nonlinear sys- 
tems, terms such as "frequency 
response" do not even have meaning, 
since liie output of a system driven by 
a sine wave is not likely to be another 
sine wave or even the same fieq^ncy. 
Fortunately, most systems in nature are 
either linear or can be made to appear 
as such by keeping the signals small 
mi linearizing around some equilibri- 
um state. So in practice, the assump- 
tion of a linear system is not a restric- 
tive one. fte's quite a lucky thing, 
too, since we have almost no mathe- 
matical tools to deal with the nonlinear 
system. 

Anyhow, when we studied Liq>hKe 
transfonns, we learned that a primary 
use is to solve dynamic systems &r &e 

state function x(t). Taking the 
form of Equation 16 gives us: 



(17) 



where P(s) is some polynomial in s 
obtained by substituting s for each 
derivative in the equation. We called 
die invert of P{s) die transfer fimction 
I^s) of the system wrote: 



Xis) = H(s) U(s) 



(18) 



Fortunately, a rich set of known 
Laplace transforms has bec^ tabulated, 

so controls engineers can usually gen- 
erate U{s) for most forcing iimctions of 
interest. Once the product H(s) U{s) 

has been written down, we can take the 
inverse transform to obtain x(t). 

hi "The Alphabet fk»m S to Z," we 
also learned how to obtain things like 
the transient response or the frequency 
response of a system from its transfer 
function. In particular, the frequency 
response of a contmuous system is 



obtained by substituting: 



S = ICO 



(19) 



into the transfer fimction. Similarly, 
we get the frequency response of a dig- 
ital sy^em hf i 



(20) 



So can we use the Fourier transform to 
do similar things? Yes and no. Of 
course -we can urate $mm ihe f Wnier 
transforms for ( Hlfe Brotial equations, 
because we can define die transform of 
a ^atfvttive. By definition: 

A.fit)\^^ j ht)e-^dt (21) 

We can integrate this equation by 
piots, using the standard formula: 

j udv = uv - j vdu 
Let: 

u = e^<^ (22) 
and: 

4v*mdt (23) 
1%aa: 

du = -icoe-'"" (24) 
and: 

Our integral becomes: 

j ht)e'^'dt = \B--^'f{t)\Z^ + 

— cc 

CO 

m j f{t)e-''"dt 



The first term, where the fimction is 
evaluated at die mfinite limits, must 

vanish. This is required by the fact that 
the total integral is finite. If the values 
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at infinity were anything but zero, this 
would not be true. So we're left with: 

oo 

7t/(0] = ^//(/)e-^^ 



(26) 



In other words, taking fte transform of 

the derivative of a function is equiva- 
lent to simply multiplying the trans- 
fonn of fbe function itself by im. 
Again, this is very much like the 
Laplace transform, where the deriva- 
tive of x(f) transforms to sX(s). 

Theoretically, we can deal with dif- 
ferential equations using Fourier trans- 
fimas bet in practice, we almost never 
do because it's simply easier to use 
Laplace transforms. In contrast to the 
rich set of known Laplace tramtfisffiaffi, 
the known Fourier transforms are lim- 
ited to those we've already discussed: 



the transforms of sine waves, delta 
fonctions, constants, and unit pulses, 
^so, Hae Courier transform has an 
awkward and difficult condition built 
in: Its integral ranges for all values of 
&<m nep^^ i[!#ltty to positive infin- 
ity. In flie teai -world. 1 can't possibly 
know die fbtoc&m over such a range, 
nffl- do I want m lOim k. Worse yet, 
the inverse transform of a function 
yields a time ftmction that extends to 
heA pc^tdve mai ^sive tee. That's 
a very unsatisfying behavior. 

We like to think that actions are 
reactions are causal. That is, a ^rst^ 
responds to some forcing function by 
moving in a certain way for all time 
after the forcing ftmction is applied. 
But the fact that the inverse transform 
extends to negative times implies that 
^ system will have resp<MHied to the 
forcing function by moving before the 
force was applied. That sort of behav- 
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ior just seems nonsensical to us. 

It was precisely the efforts to fix this 
problem that led to die definition of 
Laplace transforms. We don't have 
time or space to go into tihe details 
here, but the general idea is to multiply 
some fimction extending over all time 
by a unit step, which is zero for nega- 
tive time. That gives us a fimction 
extending only over positive time, but 
the resulting fimction, as it turns out, 
does not have a Fourier transform. 
However, we can approximate the step 
fimction by a function representmg an 
exponential decay, which does have a 
Fourier transform. We can then let the 
time constant of the exponential func- 
tion go to infinity, so that it looks like 
a unit step. The resulting formulas lead 
directly to the one-sided Laplace trans- 
form. In short, the Laplace transform is 
a direct result of tweaking the Fourier . 
transform to deal with causal systems. 
Since the Laplace transform is also 
easier to work with (more known 
transforms), we really have no incen- 
tive to use Fourier transforms for solv 
ing differential equations. 

So what good are they, then? Yo 
already know the answer — ^we di 
cussed it in the first article of thi 
series. You apply, in effect, a Fourie 
transform every time you tweak tti 
tone control knobs on your stereo. 
Remember, Fourier transforms take 
fmm the time domain to the ftequen 
domain and back. In this domain, fre- 
quency response is king. Whenever w~ 
can describe a system by its fiequen 
response, the frequency domain is the 
logical domain to operate in. An' 
whenever a signal's frequency conter 
is altered by passing it through a fre' 
quency-dependent filter, we c 
deduce the nature of the altered sig 
by transforming it to the ftequen 
domain and back again. 

The situation is not so much parall 
with that of Laplace transforms 
orthogonal to it. Recall that, wi 
Laplace transforms, we usually don 
care about the detailed time history o 
our solution; we're only interested 
the nature of the svstem. that is. its 



bility, transient response, frequency 
response, and so <hi. Wbile we can the- 
oretically solve for the detailed motion 
by taking the inverse transform, we're 
rarely motivated to do so. Laplace 
transforms are about studying the 
behavior of dynamic systems. 

In ecmtmt, Fourier transforms 
aren't very useftil for studying the 
behavior of systems, but they're great 
for studying and repnxiucing signals 
that have been filtered in some way. To 
do that, we typically need both the for- 
afii tev«ase transforms. 

Each transform has its uses. 
Conceptually, we can describe a given 
fuactkn k asy tm ^ snaMie 
domains: 

h{t)<=>H{s)<=>H{z)<=>H{(o) (27) 
Each domain has its strengths and 

FIGURE 3 

Ideal filters. 



weaknesses. The fme <ibmain is the 
one -mi ttve in, of course, so ultimately 
we'd like our results to end up here. 
Also, things Uke data sampling, gatiivg, 
and delays are best desofljad iit iMs 
domain. The s-domain is great for 
studying stability and designing con- 
trol systems. The discretC'^SNS worid 
of z is essential if we're dealing with 
digital systems. Not suiprisingly, fre- 
quency-dependent effeete are best 
described in the w-domain of Fourier 
transforms. Three ideal filters, high- 
pass, low-pass, and bandpass, are 
shown in Figure 3. The description of 
these filters is trivial in a>-space — 
ibey're Ae simple rectangles shown in 
the figure. Describing diese ideal fil- 
ters in any other domain is anything 
bttttrii^al. 

The bottom line is this: The more 
domains we have to work with, and the 
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more tools we have to work in them, 
the easier our work will be. Everything 
diat we ^ in the s, z, or frequency 
domains can be done only in the time 
domain; the Job just becomes that 
amch haider. The Fourio- ttansfonn is 
one more tool for one more domain to 
make our lives easier. We use it when 
It's useful to us. We don't use it to 
determine frequency responses of sys- 
tems because the Laplace and z trans- 
§mtm do that task better. Wh^ae Ae 
Fourier transform shines is in the 
tweaking of the frequency content of a 

THE CONVOLUTION INTEGRAL 

Until now, l\e Stadkmslj 
avoided dealing with the 
notion of a convolution inte- 
gral, mainly because we Ittvent had 

the tools to talk about it. Now that we 
have these tools, it's time to introduce 



the idea. It depends on two key con- 
(sqitS that you should by now have 
fittnly mteid: 

• The principle of superposition 

• The fflpilK. 

Siq>p08e Wm *0 pis a signal 
x(t) throng a Madt hex, and record die 

output signal y(t). We don't know yet 
what the response of the black box will 
be, but we do know h«Sw it responds to 
a unit impulse. Now suppose we tested 
it with a tmit impulse and found it to 
have a response g^t), whicft wet'd then 
use to deduce the output 

The key to the solution is die princi- 
ple of supoposition. If the system 
inside the box can be described by a 
linear system, we're safe in assuming 
diat we can treat die input signal as a 
superposition of signals. In particular, 
we can imagine it to be made up of an 



infinite number of tiny signals, each of 
which is an impulse. 

We aheady know that the response ^ 
to an impulse at f = is g(t), and that 
the response to an impulse at any other 
time r is git-r). For the impulse corre- 
sponding to time r, its height is .r(r), 
making the response of the system: 

x(t)git-v) 

The overall re^onse is represented by I 
the sum of all these litde impulse] 
responses, or: 



)<r) = /x(r)f(f-r)rfT 



(28)1 



I've extended the integral over 
time, partly because we'd like to have 
results that work for noncausal as welQ 
as causal systems and partly because 
that's the range the Fourier transform! 
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wants anyway. Equation 28 is the con- 
volution integral. If we decide to woric 
purely in the time domain, we must be 
prepared to evaluate this integral for 
each input signal. Now let's take the 
FamBmtmB^m^ of diis iat^al: 



30 



1X3 

jx(.r)g(t-r)dt 



'dt 



If we exchange the order of integra- 
tion, this becomes: 



30 



dx 




— oo 

Equation 29 tfien becomes: 



(30) 



Y(.a)) = 27iG{co) 



^JxWe-^^dr 



or. 



Y(aj) ^ 2n: G{oj)X(,w) (31) 
Except for that bothersome factor of 



2;r that most authors manage to ignore, 
the resulting transform is simply the 
product of the other two transforms. So 
we can get the response of the system 
to any input by multiplying two tran- 
forms together and then performing an 
inverse transform on the product. 

We get precisely the same kind of 
result (but without the 2;r) when we're 
dealing with Laplace transforms and z- 
transforms. If you agree that multiply- 
ing two fiinctions that are typically 
ratios of polynomials is easier than 
integrating a convolution integral, then 
you understand the value of integral 
transforms. It's precisely this: We can 
compute the responses of very com- 
plex systems to very complex signals 
amply by multiplying the individual 
ttans&r functions. 

This process of compi^^ the 
response of a system to a given signal 
goes like this: 

• Compute the transform of the signal. 

• Compute the transfer function for the 
system, multiplyii^ individual transfer 
functions if necessary. 

• Multiply the signal and the system 
transforms. 

• Take the inverse transform of the 
result. 

This process applies to all trans- 
forms. ITie two distinguishing differ- 
ences for the Fourier transform are that 
the frequency response is the transfer 
function, since we're in the frequency 
domain, and diat we tyfnlcally fmd the 
frequency response via some other 
mechanism, except when we're deal- 
ing with filters like those of Figiwe 3. 

Depending on the nature of our sys- 
tem, we may find that certain parts are 
best done in the frequency domain and 
other domains. I read of one simulation 
that required eighteen separate trans- 
formations to deal with different parts 
of the problem. In cases such as this, 
we need to be able to compute the 
transforms very efficiently. Unlike the 
Laplace transform, in which we can 
often get the desired transform from a 
table, we're almost certain to have to 
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BACK TO THE REAL WORLD 

Let's leave the world of theory 
now mi get back to practical 
implementations. The transition 
may surprise you. hronically, after all 
our efforts to #e m^A of the 

Fourier series, with its discrete coeffi- 
cients to the infinite and integral world 
of Foiiiiar tcam^^cmSt we're now 
heading back the other way. The rea- 
son is this: We can't store an infinite, 
contiouous function in ^ e^Muputra:, 
only a finite set of samples. So the inte- 
grals we worked so hard to derive must 
become simmiations again. The rou- 
tines I've given you for the Fourier 
series work just dandy for Fourier 
transforms; we need only chmige the 
interpretation of the results. 

Take a look at the forward Fourier 
transform: 

oo 

F(.a>) = ^ffit)e-^dt (32) 

— oo 

If we're going to represent the fimction 
f(t) as a table of numbers, our first 
challenge is to deal with the range, par- 
ticularly the "minus infinity" side of 
the integral. For practical reasons, 
we'd like to imagine the function 
begiiming axt = and moving forward, 
but that's not what the integral calls 
for. is it? 

I have two arguments to give that 
allow us to treat the function m begin- 
ning with time zero. All three argu- 
ments are pretty weak, and strike me 
more as arm-waving ikmx fact, but 
they're all we have. The first involves 
going back to the Fourier series, where 
die integral had a finite range: 

+71 

Cn=^f f(x)e-""dx 

The important thing to notice is that the 
t\vo end points, -n and jt. are really the 

same point: The angle is the same 
either wav, and so is the tunction of 



tiiat angle. We're not so much integrat- 
ing over a range as we are integrating 
around a circle. For that reason, it real- 
ly doesn't matter much where we 
break the circle, and changing the 
range to 0..2;r gives an equally valid 
(but difforat) r^resraitation. In fact, 
using any startiog point is equally 
valid. 

In Ae case of the Fourier transform, 
the mathematics says to integrate over 
all possible numbers, but once we 
make die commitnwnt (as we must) to 
chop the ftmction off at some point and 
make it finite at both ends, it really 
doesn't matter where we make the 
chops. 

The second argtmient is related to 
the first and is based cm the recognition 
that any finite representation of the 
function takes us right back firom a 
nonperiodic to a periodic function. The , 
idea is that we'll consider a tabular 
function to be not the whole function 
of time, but only typical of it, and we'll 
mentally extend the fimction over all 
time by tacking copies on at both ends. 
After all the trouble we went to to 
define the Fourier transform to deal 
with nonperiodic fimctions, it seems a 
bitter pill to have to recognize that 
we're back to periodic functions aga 
But that's the fact of the matter. 

This situation is very much like the 
artifice we used to arrive at the Foiuier 
transform in the first place. Recall 
we began by approximating a singl 
pulse as a square wave. Then we imag-^ 
ined stretching the waveform out with 
out changing the shape of the pulse, : 
that the periodic pulses became farthe. 
and farther apart. At the limit, w" 
imagined they'd moved all the way o" 
to infinity, where they'd be no troubl 
In fact, they're really not at infinity 
they're still hovering out there on 
horizon. The best we can hope for is to 
make the distance between them f" 
enough that the results we get 
approximate the integral transform. 
And the notion of approximation is th" 
bottom line. The instant we decide 
represent an infinitely long, continuo 
function by a finite table of poin 
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we're no longer talking about an exact 
solution but an approxinaation to it, md 
we can only hope that the a^pm^sm- 
tion is good enough. 

If you can accept this admittedly 
arm-waving argument, then we can 
represent f(t) by a tabular fimction con- 
taining N points. We'U call this func- 
tion /i^ and define ttie i^dMotlMf 
between k and t by: 



(33) 



where h is the sampling interval. 







(34) 



Take a look at that coefficient in 
front. What is its relationship to A^, the 
number of points? Well, if we accept 
the fact that we're really woricing witfi 
a periodic series again, we must also 
accept the fact that our new range is 
0..23f radians of angle, and tlmt fte 
sampling interval * tate: »^ 
up into parts. 

1^ Tbe the total time span co¥^d 
byti 



(35) 



This, of course, is the period of our 
newly periodic function. Now, where 
there's a period, there's a frequency 
given by the reciprocal of that period 
and a base angular rate: 



Ojt = iTtIT 



(36) 



Remember, we're back to periodic 
functions, so the frequency spectrum is 
no longer continuous but discrete. We 
shouldn't have a problem with that 
concept at this point, because we can 
only represent a continuous function 
by a table. As when we were defining 
the Fourier transform, we can think of 
O as being a A<b, or the difference 
between two successive fisqumcies in 
the table. 

Eliminating Tbetween Equations 35 
and 36 gives us a relationship between 



h and N: 



or: 
h = 
on 

Jl 
2n 



In 
QN 



QN 



This Is eo^elM m '. 
nMsii becomes: 



iV-l 



(37) 
(38) 

(39) 
34, 

(40) 



We've now got the first half of the 
transition from continuous to discrete 
functions. We get the second half by 
repeating the process with cv. We 
already have a big hint as to how that 
should go from the notion that Q is 
equivalent to a Aoj. We can define a set 
of discrete values of co given by: 



nQ. 



(41) 



If we substitute this defhiition into 
Equation 39, it gives us our final 
transform: 



1 

QN 



-iiiiah 



This is the official, rigorous definition 
of the discrete Fourier ffansform, and 
is i^aM for both real and eoa^lt^ ftee^ 
tions^i-. 

A word of warning: You will not 
find tilts equation in pted^y tias fetm 
in any textbook, at least none of the 
books I've researched. That's because 
people tratd to play a little fast and 
loose with minor details of the defini- 
tions. 1 trankly admit that what was 
most difficult for me in preparing this 
series was resolving seemingly mis- 
placed factors of 2.T and similar details 
between one textbook and another. 
Few authors even define the trans- 
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Itans the same way, so don't be sur- 
jnised if you find differences between 

what you see here and what you see 
elsewhere. I don't claim to be a more 
aotaii&tive source than any other ref- 
erence, but I have made considerable 
effort to keep things consistent. 

MfKUticular, you will almost always 
find the formula for the series in 
Equation 42 with the Q missing. Not 
only can tts he confusing, it isn't even 
dimffl[isionally correct. Remember, 
unlike Ihe case with the Fourier series, 
the Fourier transform and its inverse 
do not have the same units; they diifer 
by a tactor of inverse time, which is 
precise^ ttie units of Q. Folks get 
away wMi this by normalizing the 
coefiGieDts to an angular rate of one, 
correspond!^ to a period T of Iti (m 
whatever time units you'd like). If we 
do the same thing here, we'll get (using 
t3S>: 



,v-i 



-271 ink. 



/N 



(43) 



This is the form you're tm^ mme 

likely to see. It also happens to be pre- 
cisely the same form we got for the 
coefficients of the complex Fourier 
series. What does this mean? Simply 
that if we normalize the angular rate to 
OIK, we already have the code for Ae 
Fourier transform, at least the trans- 
form for real functions. It's the func- 
tion fourfit given last month md 
repeated here in Listing 1 with a mm 
appropriate name. 

One small item that remains to be 
dealt with is the range of the index n. 
Remember that in the complex Fourier 
series, the coefficients c„ range through 
negative as well as positive indices, 
with the negative indices correspond- 
ing to the negative frequencies. In 
fiinction fourfit and its companion 
fourier, we did not need or allow the 
negative- frequency coefficients be- 
cause we limited our attention to real 
functions. For real fiinctions, the nega- 
tive-frequency coefficients are the 
complex conju^tes of their positive- 



UST1NG 1 

Discrete Fourier transform. 

II Cofflpute Fourier coefficients 
fvoid fourfit(Vector 4f , CJector tc] 

int n = f.rank(); 

ijit m; -'cirankO r 
-i6,assert{Bi;<f n/2)r 
'^''awfilfir^^* exp(coiiiplex(0, 2.0 * 

r ccfflplfix'id = complexd, 0); 
'coflfiLex Mij, J ir-c 

: - for(iirt:.i=0; i < n;; i+*-){ 

for(ifrt.^;fO;f j <. m;, j++){ 

c[jl += f[i] * wij; 
. U.J' *= \a; 





frequency counterparts, so we didn't 
have to store both halves. But what 

about when we have to deal with nega- 
tive indices? C and C-H- won't let us 
use negative miMces, so how do we 
shift things around? 

Fortunately, a simple solution exists; 
it's the same one we used to define f(t). 
We haven't yet talked about a phenom- 
enon called aliasing, but we have 
talked about Shannon's limit and die 
fact that you can't get information out 
of a function at a frequency higher than 
half the sampling finquency. Aliasing 
is the reason for that, and we'll learn 
more about it next month. Briefly, the 
phenomenon causes false spectra to 
appear at intervals of the sampling fre- 
quency. Usually, this is an unwelcome 
phenomenon and something to be 
avoided. In this case, though, it's very 
welcome because it means that we 11 
find those missing n^gative-fiBquency 
coefficients repeated again at positive 
indices above ^In. The situation paral- 
lels the switch from the range -;r..;r to 
Q..2n. By extending our inctex range 



over the full N points, we pick up the 
aliased version of the negative-fre- 
quency coefficients above the limit. 

Using this concept, and following a 
similar procedure for die transition 
from integral to sum, we get the 
inverse discrete Fourier transform: 



N-l 



(44) 



or, widi die same normalized rate: 



-2.T ink. 



k/ 

/n 



(45) 



We now have forms for both the 
Fourier transform and its inverse in 
discrete form. They are: 



N-l 




Before we close, I'd like to show 
you one final form for both, in which 
we retain the rigor but in terms of T 
instead of £2. In diis form, the trans- 
form pair is: 



fk=^TFj'^/s 



T 





(46) 



N-l 



T V ^ -2;ruii / 



Again, these equations reduce to 
more common forms when we normal 
ize the period so that T= In. I^j] 
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